{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cubic interpolation\n",
    "\n",
    "### Considering unit interval [0,1]\n",
    "Let's assume that we want to interpolate the value in between two given points. And we are not interested in linear interpolation, but would like to fit a polynomial of degree 3. \n",
    "Let's also assume that we want to interpolate for values $x \\in [0,1]$, where f(0) and f(1) are given.\n",
    "Since we are fitting curved function (a polynomial), we would like to preserve the smoothness on the borders.\n",
    "I we would consider to fit a polynomial to more than 2 points. We would anyway fit a polynomial for every pair of points. To preserve the smoothness, the Cubic Hermit Spline (the thing, which we will compose below) takes into account the derivative at these two points. So, for fitting the cubic spline, we also assume that f'(0) and f'(1) are given. Either by someone mighty or compute using the finite difference. To compute the numeric derivative using finite distance we need to also know the value of the function outside [0,1].\n",
    "\n",
    "The general equation for the polynomial of degree 3 is:\n",
    "\n",
    "$$\n",
    "    f(x) = ax^3 + bx^2 + cx + d\n",
    "$$\n",
    "\n",
    "\n",
    "So to fit f(x) we need to find the coefficients (a,b,c,d). For now we have 4 unknowns (a,b,c,d) and 2 equations (for f(0) and f(1)), which gives us underdefined system of linear equations. (Note: the unknowns are a,b,c and d and they are linear).\n",
    "\n",
    "To compensate for this problem we take into account the first derivatives in given points f'(0) and f'(1).\n",
    " \n",
    "$$\n",
    "    f'(x) = 3ax^2 + 2bx + c\n",
    "$$\n",
    "\n",
    "Then we can obtain the following system of equations (1):   \n",
    "\\begin{eqnarray}\n",
    "    f(0)  &=& d \\\\\n",
    "    f(1)  &=& a + b + c +d \\\\\n",
    "    f'(0) &=& c \\\\\n",
    "    f'(1) &=& 3a + 2b + c\n",
    "\\end{eqnarray}\n",
    "\n",
    "Solving for a, b, c, d:  \n",
    "\\begin{eqnarray}\n",
    "    d &=& f(0) \\\\\n",
    "    c &=& f'(0) \\\\\n",
    "    b &=& 3f(1)-2f'(0) - 3f(0) - f'(1) \\\\\n",
    "    a &=& -2f(1) + f'(0) + 2f(0) + f'(1)\n",
    "\\end{eqnarray}\n",
    "\n",
    "Ok, found polynomial coefficients. But how do we know the derivatives given just f(0) and f(1)? This is a point in time where we need two extra points and the knowledge of finite differences :smile:\n",
    "\n",
    "Let's consider the following example. We are given points $[-1,p_0]$, $[0,p_1]$, $[1,p_2]$ and $[2, p_3]$, where $x,f(x)]$. and we want to interpolate between $[0, p_1]$ and $[1, p_2]$.\n",
    "Then,  \n",
    "\\begin{eqnarray}\n",
    "    f(0)  &=& p_1 \\\\\n",
    "    f(1)  &=& p_2 \\\\\n",
    "    f'(0) &=& (p_2 - p_0) / 1-(-1)  \\\\\n",
    "    f'(1) &=& (p_3 - p_1) / 2  \\\\\n",
    "\\end{eqnarray}\n",
    "\n",
    "the following holds:\n",
    "\\begin{eqnarray}\n",
    "    a &=& 0.5 * (3p_1 - p_0 - 3p_2 + p_3) \\\\\n",
    "    b &=& p_0 - 2.5p_1 + 2p_2 - 0.5p_3\\\\\n",
    "    c &=& 0.5(p_2 - p_0) \\\\\n",
    "    d &=& p_1\n",
    "\\end{eqnarray}\n",
    "\n",
    "The final polynom then is:\n",
    "\n",
    "\\begin{equation}\n",
    " f(x) = 0.5 * (3p_1 - p_0 - 3p_2 + p_3) x^3 + (p_0 - 2.5p_1 + 2p_2 - 0.5p_3) x^2 + 0.5(p_2 - p_0) x + p_1 \n",
    "\\end{equation}\n",
    "\n",
    "**Warning** This only works if your x is between 0 and 1 and all the points are located at distance 1 from each other in x axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Polynomial coefficients:  4.0 -6.5 0.5 5\n"
     ]
    }
   ],
   "source": [
    "## Code for testing the cubic interpolation\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.array([-1, 0, 1, 2])\n",
    "fx = np.array([2, 5, 3,4])\n",
    "\n",
    "a = 0.5 * (3 * fx[1] - fx[0] - 3*fx[2] + fx[3]);\n",
    "b = fx[0] - 2.5 * fx[1] + 2 * fx[2] - 0.5 * fx[3];\n",
    "c = 0.5 * (fx[2]- fx[0])\n",
    "d = fx[1]\n",
    "print(\"Polynomial coefficients: \",a,b,c,d)\n",
    "\n",
    "x_t = np.arange(-1, 2,0.1);\n",
    "f_x = a* np.power(x_t,3) + b*np.power(x_t,2) + c*x_t + d;\n",
    "\n",
    "# fig=plt.figure(figsize=(16, 4), dpi= 80, facecolor='w', edgecolor='k')\n",
    "fig = plt.figure(1)\n",
    "plt.plot(x,fx, 'r*')\n",
    "plt.plot(x_t,f_x, 'b-')\n",
    "plt.title(\"Cubic interpolation image\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matrix formulation\n",
    "In the above derivation, we basically solved the system of linear equation with respect to a,b,c,d by explicit substition. In general, we can rewrite this system in matrix notation and solve using technique for solving systems of linear equation using matrices. We can rewrite the Equations (1), as:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "    0 & 0 & 0 & 1 \\\\\n",
    "    1 & 1 & 1 & 1 \\\\\n",
    "    0 & 0 & 1 & 0 \\\\\n",
    "    3 & 2 & 1 & 0\n",
    "\\end{pmatrix} \n",
    "\\begin{pmatrix}\n",
    "a \\\\ b\\\\ c \\\\d \n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "f(0) \\\\ f(1) \\\\ f'(0) \\\\ f'(1)\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Afterwards solving for (a,b,c,d) should give the result. Let's check it with a code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Checking matrix formulation\n",
      "Polynomial coefficients:  4.0 -6.5 0.5 5.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XecVOX1x/HPoQpYKVFBBAsmKmo0G2KLsSU/rMRCxBIxaggmajQae6IxgSj2goVYsAIRGxpsMaixZ62AWJGOskhEI305vz/OrKzLllmm3Cnf9+s1r52duXvvuTO7Z595nueex9wdEREpHS2SDkBERLJLiV1EpMQosYuIlBgldhGREqPELiJSYpTYRURKjBK7NMjMppnZvg0890Mze28N9nmemd2SeXS5Z2YXmdndGfz8ZDPbM4sh1ez3MTMbmO39SulolXQAkltmdhTwO+A7wJfAm8AQd38+k/26+7+Bb6/Bzw1Nd1szuwjY0t2Pae5x8s3MRgKz3P2CmsfcfdtcHMvd98vFfqV0qMVewszsd8DVwFBgQ2BT4AagX5Jx5YuZqeEi5cnddSvBG7Ae8D+gfyPbjAT+Uuv7PYlWZ83304BzgXeA/wK3A2s1sG134AGgCvgMuL6BY14E3J263xNwYCAwA5gPnJ96ri+wDFieOo+3ap3XrcBcYDbwF6Bl6rnjgBeAq1Ix/KXWY9cDC4F3gX1qxdMVGAcsAD4EfllfrKnv7wM+Se3nOWDb1OODUnEuS8X6SK3Xb9/U/bbEP9k5qdvVQNvaryVwBjAvdW6/aOR9ewY4sZ5z/hyYCuyaenxman8Da/3sAcAbwBep5y+qs+9jgemp1+8Pdc6hBXAO8FHq+b8DHZP+Xddt9Zta7KVrF2At4MEM93M08H/AFsBWwAV1NzCzlsCjRELoCXQDRjfjGLsT3Tr7AH80s63d/XHik8YYd1/b3XdIbTsSWAFsCewI/AQ4sda+fkAktw2BIbUe+wjoDFwIPGBmHVPPjSaSalfgcGCome3dQJyPAb2AbwGvA/cAuPuI1P1hqVgPqudnzwd2Br4L7AD04Zuv5UbEP61uwAnAcDPboIE46voB8DbQCbg3dU7fJ16jY4DrzWzt1LZfEcl7fSLJn2RmPwUws22IT3RHAxvXiqfGKcBPgR8Rr9d/geFpxih5pMReujoB8919RYb7ud7dZ7r7AiJRHlnPNn2IP/Tfu/tX7r7Em9eH/yd3X+zubwFvEYlvNWa2IbA/cFrqOPOIluqAWpvNcffr3H2Fuy9OPTYPuNrdl7v7GOA94AAz6w7sBpydivlN4BYi8a3G3W9z9y/dfSnRmt/BzNZL8xyPBi5293nuXgX8Cfh5reeXp55f7u7jiZZ/umMYH7v77e5eDYwhPj1d7O5L3f1J4pPElqlzeMbdJ7r7Snd/GxhFJGqIf2yPuPvz7r4M+CPxiarGYOIT1axar8Hh6vIqPHpDStdnQGcza5Vhcp9Z6/50IoHX1R2YnsFxPql1fxGwdgPb9QBaA3PNrOaxFnVinFn3h4DZ7l47QdWcR1dggbt/Wee5iro7SH0qGQL0B7oAK1NPdSa6ZprSNbXvujHU+KzO69fY61DXp7XuLwZw97qPrQ1gZj8ALgF6A22ILqL7asX49evn7ovM7LNa++kBPGhmK2s9Vk18OpqdZqySB2qxl66XgKXER+eGfAW0r/X9RvVs073W/U2J/uG6ZgKb5qDlVrf06EzinDq7+/qp27r+zdkn9ZUr7Wa1/hOw6jzmAB3NbJ06z9WXpI4iBp33JbooeqYer9lvU2VS5xCJsW4M+XYvMabQ3d3XA25i1TnMBTap2dDM2hGf/GrMBPar9dqv7+5rubuSeoFRYi9R7r6Q+Cg93Mx+ambtzay1me1nZsNSm70J7G9mHc1sI+C0enb1GzPbJNUnfT7xUb+uV4mkcImZdTCztcxstyycxqdATzNrkTqnucCTwBVmtq6ZtTCzLczsR43uJfrET02df39ga2C8u88EXgT+mop5e6J/u7656+sQ/1Q+I/4Z1p22+SmweSMxjAIuMLMuZtaZeG/WeI58BtYhPqUsMbM+xD+sGmOBg8xsVzNrQ3S11P6HeBMwxMx6AKTOpSxmWBUbJfYS5u5XEHPYLyBmq8wETgYeSm1yF9GnPY1ImPUl7XtTz00lBiD/Us9xqoGDiH7cGcRg5BFZOIWaLoLPzOz11P1jiS6Empk6Y4mBvsa8Qgx6zie6Uw5395ouhiOJ1vccYqD5Qnf/Zz37uJPoPpmdOvbLdZ6/FdjGzD43s4fq/jDxulUSg5wTicHX1V7LPPg1cLGZfUn8c/l7zRPuPpkYIB1N/KP+HzE+sTS1yTVEa//J1M+/TAzcSoGxb3Y9ipQWMzuOmBq4e9KxFJvUTJrPgV7u/nHS8Uj61GIXka+Z2UGpbrsOwOXEp4tpyUYlzaXELiK19WPVwHIvYIDrY33RUVeMiEiJUYtdRKTEJHKBUufOnb1nz55JHFpEpGi99tpr8929S1PbJZLYe/bsSWVlZRKHFhEpWmY2vemt1BUjIlJylNhFREqMEruISIlRYhcRKTFpJ3Yzu83M5pnZpFqPdTSzp8zsg9TXdBcGEBGRHGlOi30ksVxZbecAT7t7L+Dp1PciIpKgtBO7uz9HrAtZWz/gjtT9O2i89reIiORBpn3sG6ZqZEOsgrNhQxua2SAzqzSzyqqqqgwPKyJSGGbMgLPPhvnzk45klawNnqYKBTVYeMbdR7h7hbtXdOnS5IVTIiJF4cor47Z4cdPb5kumif1TM9sYIPV1XuYhiYgUh88+g7/9DY4+Grp3b3r7fMk0sY8DBqbuDwQeznB/IiJF4/rrYdEiOOuspCP5puZMdxxFLJD8bTObZWYnEKud/9jMPiAW+b0kN2GKiBSWr76C666Dgw6CbbZJOppvSrsImLsf2cBT+2QpFhGRonHrrdEVc/bZSUeyOl15KiLSTMuXwxVXwO67w267JR3N6hIp2ysiUszGjIlpjsOHJx1J/dRiFxFpBne49FLYdlvYf/+ko6mfWuwiIs0wfjxMmgR33gktCrRpXKBhiYgUpksvhU03hQEDko6kYWqxi4ik6cUX4d//hquvhtatk46mYWqxi4ik6dJLoWNHOPHEpCNpnBK7iEga3nkHxo2DU06BDh2SjqZxSuwiImkYNgzatYOTT046kqYpsYuINGHmTLjnnuiC6dw56WiapsQuItKEq66K+etnnJF0JOlRYhcRacSCBTBiBBx5JPTokXQ06VFiFxFpxPDhUcmx0ErzNkaJXUSkAYsWwbXXwgEHwHbbJR1N+pTYRUQacNttsZZpIZbmbYwSu4hIPVasiNK8u+wS5XmLSVYSu5mdbmaTzWySmY0ys7WysV8RkaSMGQPTpsE554BZ0tE0T8aJ3cy6AacCFe7eG2gJFHB5HBGRxi1ZAhdcANtvDwcemHQ0zZetImCtgHZmthxoD8zJ0n5FRPLu2mujtf7UU4VbmrcxGYfs7rOBy4EZwFxgobs/WXc7MxtkZpVmVllVVZXpYUVEcmLePBgyJFrq++6bdDRrJhtdMRsA/YDNgK5ABzM7pu527j7C3SvcvaJLly6ZHlZEJCcuvDCmOV52WdKRrLlsfMjYF/jY3avcfTnwALBrFvYrIpJXkyfHVaYnnQTf+U7S0ay5bCT2GcDOZtbezAzYB5iShf2KiOTVmWfCuutGq72YZTx46u6vmNlY4HVgBfAGMCLT/YqI5NPjj8ftiiugU6eko8mMuXveD1pRUeGVlZV5P66ISH1WrIDvfjemOU6eDG3bJh1R/czsNXevaGo7rXkqImXv1lsjod9/f+Em9eYowhmaIiLZ88UX8Ic/wB57wCGHJB1NdqjFLiJlbehQqKqCxx4rvtIBDVGLXUTK1scfx+pIxx4L3/te0tFkjxK7iJStc8+Fli3jStNSosQuImXppZeiguPvfw+bbJJ0NNmlxC4iZccdTj8dNt44Enup0eCpiJSdMWPglVfg9tth7bWTjib71GIXkbKyeHEsdbfjjjFoWorUYheRsnLllTBjBowcWZy11tNRoqclIrK6t9+Giy+Gww+HvfZKOprcUWIXkbKwZAkccwxssAHceGPS0eSWumJEpCxccAFMnAj/+Ad07px0NLmlFruIlLwJE6Jv/aSTYP/9k44m95TYRaSkLVwIAwfCllsW93J3zaGuGBEpaSefDHPmwIsvQocOSUeTH1lpsZvZ+mY21szeNbMpZrZLNvYrIpKJv/8d7r47yvL26ZN0NPmTrRb7NcDj7n64mbUB2mdpvyIia2T2bBg8OBL6eeclHU1+ZZzYzWw9YA/gOAB3XwYsy3S/IiJrauVK+MUvYOlSuOsuaN066YjyKxtdMZsBVcDtZvaGmd1iZqv1ZJnZIDOrNLPKqqqqLBxWRKR+N9wATz0VC1NvtVXS0eRfNhJ7K2An4EZ33xH4Cjin7kbuPsLdK9y9okuXLlk4rIjI6qZMiYqN++0Hv/pV0tEkIxuJfRYwy91fSX0/lkj0IiJ5tWxZXF3aoUMsUF0qS901V8aJ3d0/AWaa2bdTD+0DvJPpfkVEmuvii+H112HEiKi1Xq6yNSvmFOCe1IyYqcAvsrRfEZG0vPgi/PWvcNxxcOihSUeTrKwkdnd/E6jIxr5ERJprzhwYMAA23RSuuSbpaJKnK09FpKgtXBgDpf/9Lzz7LKy7btIRJU+JXUSK1rJlcNhh8M47UbVxJ03bAJTYRaRIrVwJxx8PTz8Nd9wBP/lJ0hEVDlV3lMzNnQs/+hF88knSkUgZOfdcuOceGDq0dNcuXVNK7JK5P/8Znn8+5pqJ5MF118GwYVFf/ZzVLocUc/e8H7SiosIrKyvzflzJnpUrYU77LflgaXc+ZjNW0IpWrKAl1bRq3YKWd95Oy5bQqhVff23VCjbaKOpit1eZOFlD998P/ftDv34wdmz8fpULM3vN3Zucgag+dmmQO8ybBx98AO+/H19r3xYv/bD+H1wOHNn4vrt1iwTfq1fcau5vsYWSvjTs3/+Go4+GXXaBe+8tr6TeHErs8g1LlsDjj0cd68ceg88/X/Vcq1aReHv1gn32SSXl8dewxT+upU0bqF5WTfVRP2fFhX+muhpWrIDqar6+v3w5zJoV/xQ+/DC+Pvww1K0J160b7Lgj7LFHdN3vuGP5VeeT1b3zTrTSe/aEceOgXbukIypcSuzydTK/7774g/nf/6BTp7h677vfXdWq7tEjkvs3PPUsnNQXBg2K67jnToZezTv+woWrEv2HH8ang1dfhUcfjec7dIBdd40kv8ce8P3vw1prZeXUpUjMmRNz1du2jd/VTp2SjqiwqY+9TNVO5o88Al9+uSqZ9+8Pe+1VTxLPs08+iY/ezz0XF55MnBiPt20LO+8cSX6ffWC33ZKPVXJn4cJ4r6dOjd+FHXdMOqLkpNvHrsReZl57Da66KlrmNcn8kEPgZz+DPfcs7C6PBQti8s2zz8Yf+OuvxyBup05w8MHw05/Cj3+sj+ilZOHCeF+ffx7Gj4/3t5wpscs3TJ8O558f837XXx8OP3xVy7yQk3ljvvgiFlN48MHotlm4MAZe+/aNf1YHHAAbbJB0lLKmPv4YDjoI3nsvLkA66qikI0qeErsAMfg5dChce23Upj79dDj7bFhvvaQjy65ly6Il/+CD8NBDcc1Uq1bxKeSQQ2LQrVu3pKOUdL3wQrTUq6tjeuNeeyUdUWFIN7HrAqUStWxZVLnbYgu4/PKofPf++5HkSy2pA7RpEx/Tb7ghZt68/DKceSbMnAm/+Q107w577x2LL9Se6SOF5+67473aYIN4H5XUm0+JvcS4x0Ub22wDp50WRZFefx1GjozkVg5atIAf/CBqc7/7bkyTu/DCSPgnnhgXSR12GDzwQCx2LIVh5Uq44AL4+c9jQPzll8tzvdJsUGIvIS++GH8Q/fvHAOJjj8GTT8aUxXK29daR2N97L6ZRDh4cH/UPOww23DCS/YQJkVgkGYsWwRFHwJAh8X48/jh07Jh0VMUra4ndzFqa2Rtm9mi29inpWbo0Fu3dbTeYNg1uuQXefDMGEct1zcf6mMUc+Kuvjtb7E0/EbJoxY+Kj/6abwllnwaRJSUdaXmpqyN1/f3QbjhgRXWslKU8F87LZYv8tMCWL+5M0zJ4dvycjRsTK7B98ACecoEutm9KqVZR5vfNO+PRTGDUq5kdfdRVst110YV19dZRUkNx54w3o0wemTIlB7zPOKPHGSJ4K5mUlsZvZJsABwC3Z2J+k5/nn4Xvfixbm/fdHtbsOHZKOqvi0bx+Dy488Elc4Xntt9NOffjp07RpT7saOjYu6JHsefhh23z0S+QsvxKenktWuXZzojTdGn9+NN8b3ObroIlst9quBs4AGeynNbJCZVZpZZVXd4iDSLO7xe7HXXrEM2CuvaPHebOnSBU45BSor4x/mmWfG4HP//rHq/eDB8NJL8R7ImpkzB445JqYz9u4d4x477JB0VDk2dWpMxK+pcNe+fVQz+/jjnBwu48RuZgcC89z9tca2c/cR7l7h7hVdunTJ9LBla+lS+OUv4de/jq6EV1+FbbdNOqrStO22cMklMGNGDEIfeCDcdVfUrdlqK7joohiQlfQsWwaXXQbf/naUsjj/fHjmmZilVPI23jhaYUuWRKGjJUvi+xydfDZa7LsBB5vZNGA0sLeZ3Z2F/UodNf3pt94a08IeeSSuIpXcatky5sjfdVeMed1+ewy0XnwxfOc7UFEBV14Z74/U74knYPvtY3B6zz1h8mT4y1/KrPzDp5/GR76XX46vORxAzeqVp2a2J3Cmux/Y2Ha68rT5nn8+ygB89VUM+B1ySNIRyZw5MaPm3nuj68Ys/vEedVRMpdR0vehp+N3vYmB0yy3jorn99086quKlK09LhHtcTVm7P11JvTB07RoDrP/5T1zVe9FFkewHDYpP2P36ReL/8sukI82/RYvi2oFttolurL/+NcYslNTzQ7ViCpg7nHoqXH99FLS6+251vRQ695jCd++9MYVyzpyYk7333jHr46CDYJNNko4yd1aujHo9Z5wRhecGDIh+9VI+53xSEbAScO65MXj3u9/FH0cLfb4qKtXV0YU2blxM7fvoo3h8p50iyR98cFwVXArztt9/P8Yg7rorEnrv3rHg9J57Jh1ZaVFiL3LDhkUVxsGDoyumFP74y5l71K0ZNy5uNVMmu3ePVvzBB8diEsU0mLhgQSyheMcdMR7YogXsuy8MHBj1/bX4SfYpsRexv/0t+mkHDIjuF11FWnrmzYN//COS/JNPRp90q1bRmt9tt7jtumvMkisky5dHHZc77ohZWcuWxbTQgQNj0FilkXNLib1I3XdfFEPq2zdmEpRszQz52uLFUYTsuefiCsz//GdV1cnNNosEX5Pst902zX/0c+dGy2DMmIzmSldXxzq0EydGfKNHx+LjXbpEIj/22CjFoE+U+aHEXoRqLoLp0yfu11ykJuVl2bK42vXFFyPRv/BCTIGGmBnVuzf07BmLi/fsuer+ppvW6sr59a/h5pujOtwNNzR5TPf4XzBx4qrbpElR8rimlEKbNtFldOyx0fAo1pW3ipkSe5F56aXon+zVK67G0+wXqeEe88FfeCGS/bvvxgDljBnRoq5tQz6hJ9PoyTQ6soBqWlJNS1a2bE310QOpro6fWbmSr+8vWBBJfMGCVfvZeOMohlZz6907pi4W0xhAKVJiLyITJ8aFLZ06xSyKDTdMOiIpBitWxHTK6dOjXPP06TDtnUVM+9dUps1rx0Jfl5aspGW71rTcYD1atG5Jy5Z8fWvRIr6us04k7t69VyXyTp2SPjupT7qJXePWCZs6NWq+tGsXCzMrqUu6WrWK7pdNN4Uf/rDm0fZw0vCo49y2TfTrHJded4yUDs2MTtCcOdH9smxZJPWePZOOSEpCHmuSSGFSiz0hCxbA//1fzDB4+unovxTJigceWHV/+PDk4pDEKLEnYPnyqCPy/vuxLmmfPklHJCKlRIk9AX/6UwyS3ntv1BAREckm9bHn2YQJMHQoHH88HHlk0tGISClSYs+j+fNjSbBevWJdTRGRXFBXTJ64wwknRHJ/9FEtOi0iuZONNU+7m9kEM3vHzCab2W+zEVipufHGKPh0ySVRW0NEJFey0RWzAjjD3bcBdgZ+Y2a5mbw3d25collk83InToya6n37wm/1b09EcizjxO7uc9399dT9L4EpQG6Kd/75zzGd5OKLc7L7XFi8OAZJ118fRo7UYhkikntZ7WM3s57AjsAr2dwv7dqtKjEH0a9x442w1lqROQvYGWfEiuxPPKFyASKSH1lrP5rZ2sD9wGnu/kU9zw8ys0ozq6yqqmrezqdOjeLPNXVs27eHo4+OkncF7KGH4v/PGWdEPRgRkXzISmI3s9ZEUr/H3R+obxt3H+HuFe5e0aVLl+YdYOONoxD1kiXRSl+yJL7PYAGBXJs1K2bB7LRTzFsXEcmXbMyKMeBWYIq7X5l5SA0oosJG1dUxX33p0lipXqsgiUg+ZaOPfTfg58BEM3sz9dh57j4+C/tepYgKG11yCTz7LNx+O2y1VdLRiEi5yTixu/vzgFY8THnpJbjwwlhucuDApKMRkXKkyXdZtGhRjOl27w433aQFfkUkGSopkEWXXhoTdZ55BtZbL+loRKRcqcWeJdOnw7BhcMQRcXGsiEhSlNiz5Mwzo+vlssuSjkREyp0SexZMmABjx8I550T/uohIkpTYM7RiRRT26tEDfv/7pKMREdHgacZGjIjqjWPHRkkbEZGkqcWegQUL4A9/gL32gkMPTToaEZGgxJ6BP/4RPv8crrlGc9ZFpHAosa+hiROjcuNJJ8F22yUdjYjIKkrsa8AdTj01Fs8oojU/RKRMaPB0Ddx/f1xdOnw4dOyYdDQiIt+kFnszLV4cFyNttx0MGpR0NCIiq1OLvZkuuyzKB0yYAK306olIAVKLvRlmzIha64cfDnvumXQ0IiL1U2JvhrPOioHTyy9POhIRkYZla83Tvmb2npl9aGbnZGOfhea552DMmEjuPXokHY2ISMOyseZpS2A4sB+wDXCkmW2T6X4LSXV11IPp3h3OPjvpaEREGpeN4b8+wIfuPhXAzEYD/YB3srDvgnDfffDmm3DvvdC+fdLRiIg0LhtdMd2AmbW+n5V67BvMbJCZVZpZZVVVVRYOmx8rV8LQobD11rGIhohIocvb4Km7j3D3Cnev6NKlS74Om7FHH43yAeeeCy001CwiRSAbqWo2UHt5iU1SjxU9dxgyBDbbDI48MuloRETSk43E/h+gl5ltZmZtgAHAuCzsN3FPPw2vvhoDproYSUSKRcbpyt1XmNnJwBNAS+A2d5+ccWQFYMgQ6NoVjjsu6UhERNKXlXaou48HxmdjX4XixRej0NeVV0LbtklHIyKSPg0HNmDIEOjcWYW+RKT4KLHX4403YPx4OO006NAh6WhERJpHib0eQ4fCuuvCb36TdCQiIs2nxF7HlCmxkMbJJ8cKSSIixUaJvY5LLoF27aIbRkSkGCmx1/Lxx3DPPTFgWkQXx4qIfIMSey3DhkHLlrH0nYhIsVJiT5kzB267LS5G6rZaCTMRkeKhxJ5yxRVRd1311kWk2CmxA/Pnw003RaGvzTdPOhoRkcwosQPXXAOLFkVpXhGRYlf2if2LL+C66+DQQ2GbklrQT0TKVdkn9htugIUL4bzzko5ERCQ7yjqxL1oU1Rv79oXvfS/paEREsqOsE/vdd0NVlVrrIlJayjaxu8dMmO23h913TzoaEZHsySixm9llZvaumb1tZg+aWdGUzaqsjPK8gweDWdLRiIhkT6Yt9qeA3u6+PfA+UDQTBm++OWqtH3100pGIiGRXRond3Z909xWpb18GNsk8pNxbuBBGjYoLktZdN+loRESyK5t97McDjzX0pJkNMrNKM6usqqrK4mGb7+67Y0bM4MGJhiEikhPm7o1vYPZPYKN6njrf3R9ObXM+UAEc6k3tEKioqPDKyso1CDdz7jFg2rZt9LOLiBQLM3vN3Sua2q5VUxu4+75NHOg44EBgn3SSetJeegkmTYK//S3pSEREcqPJxN4YM+sLnAX8yN0XZSek3Lr5ZlhnHRgwIOlIRERyI9M+9uuBdYCnzOxNM7spCzHlzIIFMGYMHHMMrL120tGIiORGRi12d98yW4Hkw513wtKl8KtfJR2JiEjulM2VpzVXmu68M+ywQ9LRiIjkTkYt9mLy3HPw3nswcmTSkYiI5FbZtNhvvhnWXx9+9rOkIxERya2ySOxVVTB2LBx7LLRrl3Q0IiK5VRaJfeRIWL5cg6YiUh5KPrGvXBndMD/8oZa+E5HyUPKJ/V//go8+Ul0YESkfJZ/Yb7oJOnWCww5LOhIRkfwo6cT+ySfw8MNw3HFR9EtEpByUdGK/7TZYsQIGDUo6EhGR/CnZxF5dDSNGwN57w1ZbJR2NiEj+lGxif/JJmD5dg6YiUn5KNrHfdBN861vQr1/SkYiI5FdJJvZZs+DRR+H446FNm6SjERHJr5JM7HfcERcm/fKXSUciIpJ/JZnYR42C3XeHzTdPOhIRkfzLSmI3szPMzM2sczb2l4lJk2DyZC19JyLlK+PEbmbdgZ8AMzIPJ3OjR0OLFtC/f9KRiIgkIxst9quIBa09C/vKiHt0w+yzT8yIEREpRxkldjPrB8x297fS2HaQmVWaWWVVVVUmh21QZSVMnapuGBEpb00ujWdm/wQ2quep84HziG6YJrn7CGAEQEVFRU5a96NHQ+vWcOihudi7iEhxaDKxu/u+9T1uZtsBmwFvmRnAJsDrZtbH3T/JapRpWLkSxoyB/faLJfBERMrVGi9m7e4Tga97ss1sGlDh7vOzEFezPf88zJ4Nl12WxNFFRApHycxjHz0a2reHgw9OOhIRkWStcYu9Lnfvma19NdeKFXDffXDQQdChQ1JRiIgUhpJosT/9NMyfr9kwIiJQIol99GhYd13o2zfpSEREklf0iX3pUnjggZjiuNZaSUcjIpK8ok/sjz0GX3yhbhgRkRpFn9hHj4bOnWMJPBERKfLE/tVX8MgjUfCrdeukoxERKQwSU8etAAAGEElEQVRFndgfeQQWLVI3jIhIbUWd2EeNgq5dY1ENEREJRZvY//vfGDg94oiovy4iIqFoU+JDD8Hy5XDkkUlHIiJSWIo2sY8aFWuaVlQkHYmISGEpysQ+b16UERgwAKJisIiI1CjKxD52bNRfVzeMiMjqijKxjxoF224LvXsnHYmISOEpusQ+c2YsqqG56yIi9Su6xP73v8dXJXYRkfplnNjN7BQze9fMJpvZsGwE1ZjRo2MmzJZb5vpIIiLFKaMVlMxsL6AfsIO7LzWzbzX1M5n44AOorITLL8/lUUREilumLfaTgEvcfSmAu8/LPKSGjRkTX3/2s1weRUSkuGWa2LcCfmhmr5jZs2b2/YY2NLNBZlZpZpVVVVVrdLBu3eD446F79zUNV0Sk9Jm7N76B2T+Bjep56nxgCDABOBX4PjAG2Nyb2GlFRYVXVlauUcAiIuXKzF5z9yavt2+yj93d923kICcBD6QS+atmthLoDKxZk1xERDKWaVfMQ8BeAGa2FdAGmJ9pUCIisuYymhUD3AbcZmaTgGXAwKa6YUREJLcySuzuvgw4JkuxiIhIFhTdlaciItI4JXYRkRKjxC4iUmKU2EVESkyTFyjl5KBmVcD0NfzxzpTOlEqdS+EplfMAnUuhyuRcerh7l6Y2SiSxZ8LMKtO58qoY6FwKT6mcB+hcClU+zkVdMSIiJUaJXUSkxBRjYh+RdABZpHMpPKVyHqBzKVQ5P5ei62MXEZHGFWOLXUREGqHELiJSYgo+sZtZ/9RC2SvNrMEpQmbW18zeM7MPzeycfMaYLjPraGZPmdkHqa8bNLBdtZm9mbqNy3ecDWnqNTaztmY2JvX8K2bWM/9RpieNcznOzKpqvQ8nJhFnU8zsNjObl6qwWt/zZmbXps7zbTPbKd8xpiuNc9nTzBbWek/+mO8Y02Fm3c1sgpm9k8pdv61nm9y+L+5e0Ddga+DbwDNARQPbtAQ+AjYnasK/BWyTdOz1xDkMOCd1/xzg0ga2+1/Ssa7Jawz8GrgpdX8AMCbpuDM4l+OA65OONY1z2QPYCZjUwPP7A48BBuwMvJJ0zBmcy57Ao0nHmcZ5bAzslLq/DvB+Pb9fOX1fCr7F7u5T3P29JjbrA3zo7lM9SgmPBvrlPrpm6wfckbp/B/DTBGNprnRe49rnNxbYx8wsjzGmq1h+X5rk7s8BCxrZpB9wp4eXgfXNbOP8RNc8aZxLUXD3ue7+eur+l8AUoFudzXL6vhR8Yk9TN2Bmre9nsfoLWQg2dPe5qfufABs2sN1aqYW/XzazQkn+6bzGX2/j7iuAhUCnvETXPOn+vhyW+pg81syKdQn1YvnbSNcuZvaWmT1mZtsmHUxTUt2ROwKv1Hkqp+9LpisoZUVjC2a7+8P5jicTTSz+/TV3dzNraK5pD3efbWabA/8ys4nu/lG2Y5VGPQKMcvelZvYr4pPI3gnHVO5eJ/42/mdm+xNLc/ZKOKYGmdnawP3Aae7+RT6PXRCJ3RtZMDtNs4HaLapNUo/lXWPnYmafmtnG7j439bFrXgP7mJ36OtXMniH+4yed2NN5jWu2mWVmrYD1gM/yE16zNHku7l477luI8ZFiVDB/G5mqnRzdfbyZ3WBmnd294IqDmVlrIqnf4+4P1LNJTt+XUumK+Q/Qy8w2M7M2xMBdwcwmqWUcMDB1fyCw2qcRM9vAzNqm7ncGdgPeyVuEDUvnNa59focD//LUSFGBafJc6vR3Hkz0kxajccCxqVkYOwMLa3UHFhUz26hmzMbM+hD5q+AaDqkYbwWmuPuVDWyW2/cl6RHkNEaYDyH6n5YCnwJPpB7vCoyvM8r8PtGyPT/puBs4l07A08AHwD+BjqnHK4BbUvd3BSYSMzUmAickHXdjrzFwMXBw6v5awH3Ah8CrwOZJx5zBufwVmJx6HyYA30k65gbOYxQwF1ie+js5ARgMDE49b8Dw1HlOpIGZZYVwS+NcTq71nrwM7Jp0zA2cx+6AA28Db6Zu++fzfVFJARGRElMqXTEiIpKixC4iUmKU2EVESowSu4hIiVFiFxEpMUrsIiIlRoldRKTE/D8eykKCglC78wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## considering x, fx from previous cell\n",
    "from numpy import linalg as LA\n",
    "print(\"Checking matrix formulation\")\n",
    "F = np.matrix([[0,0,0,1],[1,1,1,1], [0,0,1,0], [3,2,1,0]])\n",
    "b = np.array([fx[1], fx[2], 0.5 * (fx[2] - fx[0]), 0.5 * (fx[3]-fx[1])])\n",
    "abcd = np.linalg.solve(F,b)\n",
    "print(\"Polynomial coefficients: \",abcd[0], abcd[1], abcd[2], abcd[3])\n",
    "\n",
    "f_x = abcd[0]* np.power(x_t,3) + abcd[1]*np.power(x_t,2) + abcd[2]*x_t + abcd[3];\n",
    "\n",
    "# fig=plt.figure(figsize=(16, 4), dpi= 80, facecolor='w', edgecolor='k')\n",
    "fig = plt.figure(1)\n",
    "plt.plot(x,fx, 'r*')\n",
    "plt.plot(x_t,f_x, 'b-')\n",
    "plt.title(\"Cubic interpolation image\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last equation derived above is also confirmed by Wikipedia (https://en.wikipedia.org/wiki/Cubic_Hermite_spline) on the unit interval. If we rearrange, the $f(x)$ with respect to $p_0,p_1, p_2,p_3$ rather than $x$, we get exactly the formula from Wikipedia (Checked!).\n",
    "\n",
    "**Side note** (for curious)\n",
    "After the rearrangement the polynomial will have the following form:\n",
    "\n",
    "$$\n",
    "    f(x) = (2x^3 - 3x^2 + 1)f(0) + (x^3-2x^2+x)f'(0) + (-2x^3+3x^2)f(1) + (x^3-x^2)f'(1)\n",
    "$$\n",
    "\n",
    "The \"coeffients\" in this case are called Hermite basis functions and have further properties (not relevant here).\n",
    "\n",
    "$$\n",
    "    f(x) = h_{00}f(0) + h_{01}f'(0) + h_{10}f(1) + h_{11}f'(1)\n",
    "$$\n",
    "\n",
    "### Arbitrary interval\n",
    "\n",
    "Ok. What if the x coordinates are not 0 and 1? Since, all the nice equations above hold only if we compute f'(0) and f'(1).\n",
    "I guess, there 2 cases to consider:\n",
    "    1) Your points are at 1 distance but shifted somewhere. Then the solution would be to shift the result.\n",
    "    2) The points are not equidistant or the distance is not 1. Then one would need to take distance into account for computing derivatives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bicubic interpolation\n",
    "\n",
    "[See https://en.wikipedia.org/wiki/Bicubic_interpolation for explicit computation of all terms.]\n",
    "\n",
    "Analogously to the cubic interpolation, we would like to fit a 3 degree polynomial, but now for the 2D data.\n",
    "The general equation for the 2D polynomial of 3 order is:\n",
    "\n",
    "$$\n",
    "    f(x,y) = \\sum_{i=0}^3\\sum_{j=0}^3 a_{ij}x^iy^j\n",
    "$$\n",
    "\n",
    "Rewritting this by opening the sum will give us $16$ coeffients.\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    f(x,y) =  a_{00}x^0y^0 + a_{01}x^0y^1 + a_{02}x^0y^2 + a_{03}x^0y^3 + \\\\\n",
    "               a_{10}x^1y^0 + a_{11}x^1y^1 + a_{12}x^1y^2 + a_{13}x^1y^3 + \\\\\n",
    "               a_{20}x^2y^0 + a_{21}x^2y^1 + a_{22}x^2y^2 + a_{23}x^2y^3 + \\\\\n",
    "               a_{30}x^3y^0 + a_{31}x^3y^1 + a_{32}x^3y^2 + a_{33}x^3y^3 + \\\\            \n",
    "\\end{eqnarray}\n",
    "\n",
    "As in 1D case, we assume unit square $[0,1] x [0,1]$ and to fit a polynomial, we assume the function value on the boundaries to be given f(0,0), f(0,1), f(1,0), f(1,1). Then, we estimate first 4 equations, by evaluating f(x,y) at the boundaries.\n",
    "As well as in 1D case to preserve the smoothness on the borders, we also assume to know the partial derivatives $f_{x}(x,y)$, $f_{y}(x,y)$ and mixed partial derivatives $f_{xy}(x,y)$. Evaluating these partial derivatives on the boundaries gives us further 12 equations.\n",
    "\n",
    "As in matrix formulation for 1D case, we can form a 16x16 matrix of parameters and respective f(.,.) terms and solve the system of linear equations with respect to $a_{ij}$ parameters. See [https://en.wikipedia.org/wiki/Bicubic_interpolation](Bicubic interpolation) for full matrix.\n",
    "\n",
    "In more favourable matrix form, the coeffients can be found as:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "    a_{00} & a_{01} & a_{02} & a_{03} \\\\ \n",
    "    a_{10} & a_{11} & a_{12} & a_{13} \\\\\n",
    "    a_{20} & a_{21} & a_{22} & a_{23} \\\\\n",
    "    a_{30} & a_{31} & a_{32} & a_{33} \n",
    "\\end{pmatrix} \n",
    "=\n",
    "\\begin{pmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\ \n",
    "    0 & 0 & 1 & 0 \\\\\n",
    "    -3 & 3 & -2 & -1 \\\\\n",
    "    2 & -2 & 1 & 1 \n",
    "\\end{pmatrix}  \n",
    "\\begin{pmatrix}\n",
    "    f(0,0) & f(0,1) & f_y(0,0) & f_y(0,1) \\\\ \n",
    "    f(1,0) & f(1,1) & f_y(1,0) & f_y(1,1) \\\\\n",
    "    f_x(0,0)& f_x(0,1) & f_{xy}(0,0) & f_{xy}(0,1) \\\\\n",
    "    f_x(1,0) & f_x(1,1) & f_{xy}(1,0) & f_{xy}(1,1) \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "    1 & 0 & -3 & 2 \\\\ \n",
    "    0 & 0 & 3 & -2 \\\\\n",
    "    0 & 1 & -2 & 1 \\\\\n",
    "    0 & 0 & -1 & 1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "The final interpolated value can be obtained as:\n",
    "\n",
    "$$\n",
    " f(x,y) = \n",
    " \\begin{pmatrix}\n",
    "  1 & x & x^2 & x^3\n",
    " \\end{pmatrix}\n",
    " \\begin{pmatrix}\n",
    "    a_{00} & a_{01} & a_{02} & a_{03} \\\\ \n",
    "    a_{10} & a_{11} & a_{12} & a_{13} \\\\\n",
    "    a_{20} & a_{21} & a_{22} & a_{23} \\\\\n",
    "    a_{30} & a_{31} & a_{32} & a_{33} \n",
    "\\end{pmatrix} \n",
    "\\begin{pmatrix}\n",
    "  1 \\\\ y \\\\ y^2 \\\\ y^3\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "## Finding derivatives\n",
    "\n",
    "To solve the anove specified system, we either need to know all the derivatives beforehand, or we can approximate them using the [finite differences](https://en.wikipedia.org/wiki/Finite_difference).\n",
    "The partial derivative then can be computed as:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    f_x(x,y) &=& 0.5 * (f(x+h, y) - f(x-h,y)) \\\\\n",
    "    f_y(x,y) &=& 0.5 * (f(x, y+k) - f(x, y-k)) \\\\\n",
    "    f_{xy}(x,y) &=& \\frac{f(x+h,y+k) - f(x+h, y-k) - f(x-h,y+k) + f(x-h, y-k)}{4hk}\n",
    "\\end{eqnarray}\n",
    "\n",
    "For the case of bicubic interpolation, all our values are equidistance with distance 1 (since we work with pixels). So, h and k are equal to 1.\n",
    "Considering the image patch like in the image\n",
    "![image](data/bicubic_int.png)\n",
    "If we want to find the partial derivatives for the point $c_{11}$, they will involve the following pixels:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    f_x(x,y) &=& 0.5 * (c_{12} - c_{10}) \\\\\n",
    "    f_y(x,y) &=& 0.5 * (c_{21} - c_{01}) \\\\\n",
    "    f_{xy}(x,y) &=& \\frac{c_{22}- c_{02} - c_{20} + c_{00}}{4}\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bicubicInterpolation(image, scale):\n",
    "    orig_rows = image.shape[0];\n",
    "    orig_cols = image.shape[1];\n",
    "    img_sc = np.zeros((orig_rows*scale, orig_cols*scale));\n",
    "#     append 2 columns for each border\n",
    "    B = np.matrix([[1.0, 0, 0, 0], [0, 0, 1, 0], [-3, 3, -2, -1], [2, -2, 1, 1]])\n",
    "    C = np.matrix([[1.0, 0, -3, 2], [0, 0, 3, -2], [0, 1, -2, 1], [0, 0, -1, 1]])\n",
    "    \n",
    "    for r in range(0, img_sc.shape[0]):\n",
    "        for c in range (0, img_sc.shape[1]):   \n",
    "            orig_row = int(np.floor(r/scale))\n",
    "            orig_col = int(np.floor(c/scale))\n",
    "            patch = np.zeros((4,4));\n",
    "            ## make a patch to avoid border conditions check\n",
    "            pr = 0;\n",
    "            pc = 0;\n",
    "            for y in range(orig_row-1, orig_row + 3):\n",
    "                for x in range(orig_col - 1, orig_col +3):\n",
    "                    x = max(0, x);\n",
    "                    x = min(x, orig_cols - 1);\n",
    "                    y = max(0, y);\n",
    "                    y = min(y, orig_rows - 1)\n",
    "                    patch[pr,pc] = image[y,x]\n",
    "                    pc += 1;\n",
    "                pr += 1;\n",
    "                pc = 0;\n",
    "            \n",
    "            F = np.zeros((4,4));\n",
    "            ## where to find f values in a constructed patch\n",
    "            f00 = [1,1];\n",
    "            f01 = [1,2];\n",
    "            f10 = [2,1];\n",
    "            f11 = [2,2];\n",
    "            \n",
    "            F[0,0] = patch[f00[0], f00[1]];\n",
    "            F[0,1] = patch[f01[0], f01[1]];\n",
    "            F[1,0] = patch[f10[0], f10[1]];\n",
    "            F[1,1] = patch[f11[0], f11[1]];\n",
    "            ## y derivatives\n",
    "            F[0,2] = 0.5 * (patch[f00[0], f00[1] + 1] - patch[f00[0], f00[1] - 1])\n",
    "            F[0,3] = 0.5 * (patch[f01[0], f01[1] + 1] - patch[f01[0], f01[1] - 1])\n",
    "            F[1,2] = 0.5 * (patch[f10[0], f10[1] + 1] - patch[f10[0], f10[1] - 1])\n",
    "            F[1,3] = 0.5 * (patch[f11[0], f11[1] + 1] - patch[f11[0], f11[1] - 1])\n",
    "            ## x derivatives\n",
    "            F[2,0] = 0.5 * (patch[f00[0] +1, f00[1]] - patch[f00[0] -1, f00[1]])\n",
    "            F[2,1] = 0.5 * (patch[f01[0] +1, f01[1]] - patch[f01[0] -1, f01[1]])\n",
    "            F[3,0] = 0.5 * (patch[f10[0] +1, f10[1]] - patch[f10[0] -1, f10[1]])\n",
    "            F[3,1] = 0.5 * (patch[f11[0] +1, f11[1]] - patch[f11[0] -1, f11[1]])\n",
    "            # mixed derivatives\n",
    "            F[2,2] = 0.25* (patch[f00[0]+1, f00[1]+1] - patch[f00[0]+1, f00[1]-1] - \n",
    "                            patch[f00[0]-1, f00[1]+1] + patch[f00[0]-1, f00[1]-1])\n",
    "            F[2,3] = 0.25* (patch[f01[0]+1, f01[1]+1] - patch[f01[0]+1, f01[1]-1] - \n",
    "                            patch[f01[0]-1, f01[1]+1] + patch[f01[0]-1, f01[1]-1])\n",
    "            F[3,2] = 0.25* (patch[f10[0]+1, f10[1]+1] - patch[f10[0]+1, f10[1]-1] - \n",
    "                            patch[f10[0]-1, f10[1]+1] + patch[f10[0]-1, f10[1]-1])\n",
    "            F[3,3] = 0.25* (patch[f11[0]+1, f11[1]+1] - patch[f11[0]+1, f11[1]-1] - \n",
    "                            patch[f11[0]-1, f11[1]+1] + patch[f11[0]-1, f11[1]-1])\n",
    "            A = B.dot(F).dot(C);\n",
    "            dx = c/scale - np.floor(c/scale);\n",
    "            dy = r/scale - np.floor(r/scale);\n",
    "            Y = np.matrix([[1, dy, np.power(dy,2), np.power(dy,3)]])\n",
    "            X = np.matrix([[1], [dx], [np.power(dx,2)], [np.power(dx,3)]])\n",
    "            img_sc[r,c] = Y*A*X;\n",
    "            \n",
    "    return img_sc;\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[10  0]\n",
      " [ 0  0]]\n",
      "[[10.         5.         0.        -0.625    ]\n",
      " [ 5.         2.5        0.        -0.3125   ]\n",
      " [ 0.         0.         0.         0.       ]\n",
      " [-0.625     -0.3125     0.         0.0390625]]\n"
     ]
    }
   ],
   "source": [
    "## testing the interpolation\n",
    "m = np.matrix([[10,0], [0,0]]);\n",
    "scale = 2;\n",
    "print(m)\n",
    "img_sc= bicubicInterpolation(m, scale);\n",
    "print(img_sc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAADT5JREFUeJzt3X+o3fV9x/Hnq0nUoa3axGFI4o+hSDu3ac0yi7CJVlDXmcGUKaU1RclWqtWxQtsNHC2M2cJa6CwtQWUqpbVo57LikAzt2rLFGUNMNc42E4axUm20sWlt7HXv/XG+uuvtTT6Z53u/517zfMDhfr/nfDzv9yHJy3O/33O+71QVknQgb5l0A5LmP4NCUpNBIanJoJDUZFBIajIoJDWNFRRJ3p5kU5Lvdz+P3c+6V5Js624bx6kpaXgZ53MUST4DPF9VNyb5OHBsVX1slnV7q+qoMfqUNEHjBsUTwLlV9UyS5cA3q+q0WdYZFNICNm5Q/Liqjum2A7zw6v6MdVPANmAKuLGq7tnP860H1gO8ZfFhZx3xtl99w73NV4v3vjzpFubM1FGHTbqFOfGOlc9NuoU58/D2fT+qquNa6xa3FiT5F+D4WR76y+k7VVVJ9pc6J1bV00l+Dbg/yXer6r9mLqqqDcAGgCOXrqpf//3rW+0tOEv/ddekW5gzu39v5aRbmBObP/OlSbcwZxYt3/nfB7OuGRRV9Z79PZbkh0mWT/vV49n9PMfT3c8nk3wTOBP4paCQND+Ne3p0I3Blt30l8I8zFyQ5Nsnh3fYy4Bxgx5h1JQ1o3KC4EbggyfeB93T7JFmd5OZuzTuALUkeAR5gdIzCoJAWkOavHgdSVbuB82e5fwtwdbf9b8BvjFNH0mT5yUxJTQaFpCaDQlKTQSGpyaCQ1GRQSGoyKCQ1GRSSmgwKSU0GhaQmg0JSk0EhqcmgkNRkUEhqMigkNRkUkpoMCklNBoWkpl6CIsmFSZ5IsrObGDbz8cOT3Nk9/mCSk/qoK2kYYwdFkkXAF4CLgHcCVyR554xlVzEaDnQK8Dng0+PWlTScPt5RrAF2VtWTVfUy8FVg7Yw1a4Hbuu27gPO7yWKSFoA+gmIF8NS0/V3dfbOuqaopYA+wtIfakgYwrw5mJlmfZEuSLVP7fjrpdiR1+giKp4FV0/ZXdvfNuibJYuBoYPfMJ6qqDVW1uqpWLz78yB5ak9SHPoLiIeDUJCcnOQy4nNGowemmjx68FLi/xhmjLmlQY00Kg9ExhyTXAPcBi4Bbq+qxJJ8CtlTVRuAW4I4kO4HnGYWJpAVi7KAAqKp7gXtn3HfDtO2fA5f1UUvS8ObVwUxJ85NBIanJoJDUZFBIajIoJDUZFJKaDApJTQaFpCaDQlKTQSGpyaCQ1GRQSGoyKCQ1GRSSmgwKSU0GhaQmg0JSk0EhqcmgkNQ01OzRdUmeS7Ktu13dR11Jwxj74rrTZo9ewGhK2ENJNlbVjhlL76yqa8atJ2l4fVyF+7XZowBJXp09OjMo/l9OXvFDbv/rv+2hvfnl2hPPmXQLc+alZSdMuoU5cfrm9026hTn0yYNaNdTsUYA/SrI9yV1JVs3y+OtGCr7w/P/00JqkPgx1MPOfgJOq6jeBTfzfZPPXmT5S8Ni3e5xVmi8GmT1aVbural+3ezNwVg91JQ1kkNmjSZZP270EeLyHupIGMtTs0Y8kuQSYYjR7dN24dSUNZ6jZo58APtFHLUnD84ihpCaDQlKTQSGpyaCQ1GRQSGoyKCQ1GRSSmgwKSU0GhaQmg0JSk0EhqcmgkNRkUEhqMigkNRkUkpoMCklNBoWkJoNCUlNfIwVvTfJskkf383iSfL4bObg9ybv6qCtpGH29o/h74MIDPH4RcGp3Ww98sae6kgbQS1BU1bcYXV17f9YCt9fIZuCYGZfwlzSPDXWM4qDGDjpSUJqf5tXBTEcKSvPTUP8am2MHJc1fQwXFRuAD3dmPs4E9VfXMQLUljamXSWFJvgKcCyxLsgv4K2AJQFV9idEUsYuBncDPgA/2UVfSMPoaKXhF4/ECPtxHLUnD84ihpCaDQlKTQSGpyaCQ1GRQSGoyKCQ1GRSSmgwKSU0GhaQmg0JSk0EhqcmgkNRkUEhqMigkNRkUkpoMCklNBoWkJoNCUtNQIwXPTbInybbudkMfdSUNo5drZjIaKXgTcPsB1ny7qt7bUz1JAxpqpKCkBayvdxQH491JHgF+AHy0qh6buSDJekZDjDnisKO59rI/HbC9YSw69aeTbmHOvHja1KRbmBMP/vaGSbcwZw52APBQQbEVOLGq9ia5GLiH0WTz16mqDcAGgLcdtaIG6k1SwyBnParqxara223fCyxJsmyI2pLGN0hQJDk+SbrtNV3d3UPUljS+oUYKXgp8KMkU8BJweTc9TNICMNRIwZsYnT6VtAD5yUxJTQaFpCaDQlKTQSGpyaCQ1GRQSGoyKCQ1GRSSmgwKSU0GhaQmg0JSk0EhqcmgkNRkUEhqMigkNRkUkpoMCklNBoWkprGDIsmqJA8k2ZHksSTXzbImST6fZGeS7UneNW5dScPp45qZU8CfV9XWJG8FHk6yqap2TFtzEaM5HqcCvwN8sfspaQEY+x1FVT1TVVu77Z8AjwMrZixbC9xeI5uBY5Ic7JAiSRPW6zGKJCcBZwIPznhoBfDUtP1d/HKYkGR9ki1JtvziF2/e0XvSQtNbUCQ5CrgbuL6qXnwjz1FVG6pqdVWtXrLkyL5akzSmXoIiyRJGIfHlqvr6LEueBlZN21/Z3SdpAejjrEeAW4DHq+qz+1m2EfhAd/bjbGBPVT0zbm1Jw+jjrMc5wPuB7ybZ1t33F8AJ8NpIwXuBi4GdwM+AD/ZQV9JAxg6KqvoOkMaaAj48bi1Jk+EnMyU1GRSSmgwKSU0GhaQmg0JSk0EhqcmgkNRkUEhqMigkNRkUkpoMCklNBoWkJoNCUpNBIanJoJDUZFBIajIoJDUZFJKahhopeG6SPUm2dbcbxq0raThDjRQE+HZVvbeHepIGNtRIQUkLWB/vKF5zgJGCAO9O8gjwA+CjVfXYLP/9emA9wAkrFnPfPXf02d68cPrm9026hTlzJG/OMZDrzviDSbcwhzYc1KqhRgpuBU6sqt8C/g64Z7bnmD5S8Lili/pqTdKYBhkpWFUvVtXebvteYEmSZX3UljT3BhkpmOT4bh1J1nR1d49bW9IwhhopeCnwoSRTwEvA5d30MEkLwFAjBW8Cbhq3lqTJ8JOZkpoMCklNBoWkJoNCUpNBIanJoJDUZFBIajIoJDUZFJKaDApJTQaFpCaDQlKTQSGpyaCQ1GRQSGoyKCQ1GRSSmgwKSU19XFz3iCT/keSRbqTgJ2dZc3iSO5PsTPJgN/9D0gLRxzuKfcB53cyOM4ALk5w9Y81VwAtVdQrwOeDTPdSVNJA+RgrWqzM7gCXdbeYVttcCt3XbdwHnv3r5fknzX18DgBZ1l+p/FthUVTNHCq4AngKoqilgD7C0j9qS5l4vQVFVr1TVGcBKYE2S09/I8yRZn2RLki3P7X6lj9Yk9aDXsx5V9WPgAeDCGQ89DawCSLIYOJpZJoU5e1San/o463FckmO67V8BLgD+c8ayjcCV3falwP1OCpMWjj5GCi4HbkuyiFHwfK2qvpHkU8CWqtrIaDbpHUl2As8Dl/dQV9JA+hgpuB04c5b7b5i2/XPgsnFrSZoMP5kpqcmgkNRkUEhqMigkNRkUkpoMCklNBoWkJoNCUpNBIanJoJDUZFBIajIoJDUZFJKaDApJTQaFpCaDQlKTQSGpyaCQ1GRQSGoaavbouiTPJdnW3a4et66k4fRxFe5XZ4/uTbIE+E6Sf66qzTPW3VlV1/RQT9LA+rgKdwGt2aOSFrD0MYenm+nxMHAK8IWq+tiMx9cBfwM8B3wP+LOqemqW51kPrO92TwOeGLu5g7cM+NGA9Ybi61p4hnxtJ1bVca1FvQTFa082mhj2D8C1VfXotPuXAnural+SPwH+uKrO661wD5JsqarVk+6jb76uhWc+vrZBZo9W1e6q2tft3gyc1WddSXNrkNmjSZZP270EeHzcupKGM9Ts0Y8kuQSYYjR7dF0Pdfu2YdINzBFf18Iz715br8coJL05+clMSU0GhaSmQz4oklyY5IkkO5N8fNL99CXJrUmeTfJoe/XCkWRVkgeS7Oi+MnDdpHvqw8F8FWKSDuljFN0B2O8xOlOzC3gIuKKqdky0sR4k+V1Gn5i9vapOn3Q/fenOoC2vqq1J3srog35/uND/zJIEOHL6VyGA62b5KsREHOrvKNYAO6vqyap6GfgqsHbCPfWiqr7F6AzTm0pVPVNVW7vtnzA61b5isl2Nr0bm7VchDvWgWAFM/yj5Lt4Ef+kOFUlOAs4EHpxsJ/1IsijJNuBZYFNVzZvXdagHhRaoJEcBdwPXV9WLk+6nD1X1SlWdAawE1iSZN78yHupB8TSwatr+yu4+zWPd7/B3A1+uqq9Pup++7e+rEJN0qAfFQ8CpSU5OchhwObBxwj3pALqDfrcAj1fVZyfdT18O5qsQk3RIB0VVTQHXAPcxOij2tap6bLJd9SPJV4B/B05LsivJVZPuqSfnAO8Hzpt2xbSLJ91UD5YDDyTZzuh/YJuq6hsT7uk1h/TpUUkH55B+RyHp4BgUkpoMCklNBoWkJoNCUpNBIanJoJDU9L9yyhvGDjGp5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "img_name = 'data/inter_test.png';\n",
    "image = cv2.imread(img_name, 0);\n",
    "image = image.astype('float')\n",
    "# new_img = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n",
    "# image = image.astype('float')\n",
    "plt.imshow(image);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison between image interpolations\n",
    "Now I want to see the difference between the nearest neighbour, bilinear and cubic interpolation.\n",
    "\n",
    "### Nearest neighbour interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nnInterpolation(image, scale):\n",
    "    orig_rows = image.shape[0];\n",
    "    orig_cols = image.shape[1];\n",
    "    img_sc = np.zeros((orig_rows*scale, orig_cols*scale));\n",
    "    for r in range(0,img_sc.shape[0]):\n",
    "        for c in range(0, img_sc.shape[1]):\n",
    "            orig_pix_r = int(np.floor(r/scale));\n",
    "            orig_pix_c = int(np.floor(c/scale));\n",
    "            img_sc[r,c] = image[orig_pix_r, orig_pix_c];\n",
    "    return img_sc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bilinear interpolation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bilinearInterpolation(image, scale):\n",
    "    orig_rows = image.shape[0];\n",
    "    orig_cols = image.shape[1];\n",
    "    img_sc = np.zeros((orig_rows*scale, orig_cols*scale));\n",
    "    \n",
    "    for r in range(0, img_sc.shape[0]):\n",
    "        for c in range(0, img_sc.shape[1]):\n",
    "            orig_pix_r = r / scale;\n",
    "            orig_pix_c = c / scale;\n",
    "\n",
    "            ## establishing neighboring pixels coordinates\n",
    "            p_row_up = int(max(0, np.floor(orig_pix_r)));\n",
    "            p_row_down = int(min(orig_rows-1, np.ceil(orig_pix_r)));\n",
    "            p_col_left = int(max(0, np.floor(orig_pix_c)));\n",
    "            p_col_right = int(min(orig_cols-1, np.ceil(orig_pix_c)));\n",
    "            a00 = image[p_row_up, p_col_left];\n",
    "            a01 = image[p_row_up, p_col_right];\n",
    "            a10 = image[p_row_down, p_col_left];\n",
    "            a11 = image[p_row_down, p_col_right];\n",
    "            dx = orig_pix_r - p_row_up;\n",
    "            dy = orig_pix_c - p_col_left;\n",
    "            img_sc[r,c] = a00 * (1.- dx) * (1. - dy) + \\\n",
    "                            a01 * (1.- dx) * dy + a10 * dx * (1.- dy) + \\\n",
    "                            a11 * dx * dy;\n",
    "    return img_sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABRoAAAFPCAYAAAAiKXJ8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xl8FPX9P/DXzO5mNwES5Aw3EZBDFIQKoiigKEWlokjFi6OIIqBS1P7AKgoe4FWhiHgCXlgLii2lUPHCavFCqHhWLKBFDgFJIOTY3ZnfH/lm9/357M4w2dzh9Xw8eDAzn8/OTjab924m+3qPYdu2DSIiIiIiIiIiIqJyMKv7AIiIiIiIiIiIiKj244lGIiIiIiIiIiIiKjeeaCQiIiIiIiIiIqJy44lGIiIiIiIiIiIiKjeeaCQiIiIiIiIiIqJy44lGIiIiIiIiIiIiKjeeaCQiIiIiIiIiIqJy44lGIiIiIiIiIiIiKjeeaCQiIiIiIiIiIqJy44lGqlB33XUXDMNI6bZLly6FYRjYvn17xR6UsH37dhiGgaVLl7rOe+edd2AYBt55551KOxYiomOd15rsdtuHHnqo4g+MiGo0wzBw1113xdaTvYccOHAgBg4cWOXHRkRU1fSaWNHat2+PCy+88Kjz+Ds0leKJRgIAfPHFF7jqqqvQqlUrBINBtGzZEldeeSW++OKL6j40IqIqVfoLaygUws6dOxPGBw4ciO7du1fDkVW8ZcuWYd68edV9GEREsdor/zVr1gyDBg3CmjVrqvvwiIiqFGsi1Wb+6j4Aqn6vvvoqLr/8cjRq1Ajjx49HTk4Otm/fjmeeeQYrVqzAn/70J1x88cWe9nX77bdj+vTpKR3H1VdfjVGjRiEYDKZ0+4p01llnoaCgAGlpadV9KERUTYqKijB37lwsWLCgug+l0ixbtgyff/45pk6dWi33365dOxQUFCAQCFTL/RNRzTN79mzk5OTAtm3s2bMHS5cuxfnnn49Vq1bFPlFTUFAAv9/915jXX3+9Kg6XiKhSVVRNrAr8HZpKVf+zkarVd999h6uvvhrHH3883n33XTRt2jQ2dtNNN+HMM8/E1Vdfjc8++wzHH3+8437y8/NRr149+P3+lIucz+eDz+dL6bYVzTRNhEKh6j4MIqpGPXv2xFNPPYUZM2agZcuW1X04sTpbl5R+cpSIqNTQoUPxi1/8IrY+fvx4NG/eHC+99FLsl2ovdaOm/6JbF2s6EVW8iqqJVYG/Q1MpRqePcQ8++CCOHDmCJ598UjnJCABNmjTBE088gfz8fDzwwAOx7aV9GL/88ktcccUVOO6449C/f39lTCooKMCNN96IJk2aoEGDBvjVr36FnTt3euqvU9oP4r333kOfPn0QCoVw/PHH47nnnlPu48CBA7jllltw0kknoX79+sjMzMTQoUPx73//O6XHJVl/idK45GeffYYBAwYgIyMDHTt2xIoVKwAA69evR9++fZGeno7OnTvjjTfeUPa5Y8cOTJo0CZ07d0Z6ejoaN26MkSNHJu1JWXof6enpaN26Ne655x4sWbIkaQ/LNWvW4Mwzz0S9evXQoEEDXHDBBYy8E1WA2267DdFoFHPnzvU0/4UXXkDv3r2Rnp6ORo0aYdSoUfjhhx+UOf/85z8xcuRItG3bFsFgEG3atMFvf/tbFBQUKPPGjh2L+vXr47vvvsP555+PBg0a4Morr4yNf/jhh/jlL3+JrKwsZGRkYMCAAXj//feVfRw6dAhTp05F+/btEQwG0axZM5x77rn49NNPAZTUtNWrV2PHjh2xSE779u1dv0bDMDBlyhS89tpr6N69O4LBIE488USsXbs2Ye7OnTvxm9/8Bs2bN4/NW7x4sTLHqUfj8uXL0a1bN4RCIXTv3h0rV67E2LFjHY/vySefRIcOHRAMBnHqqafi448/dv06iKj2aNiwIdLT05U/ZHvpR6b3aCx9b/fnP/8Z9957L1q3bo1QKIRzzjkHW7duTbi9lzrr9b1d6Xvc9evXY9KkSWjWrBlat25d5seCiMhrTdy5cyfGjx+Pli1bIhgMIicnB9dffz2Ki4sBOF9bwe26Ca+//jp69uyJUCiEbt264dVXX1XGnXo0fvjhhzj//PNx3HHHoV69ejj55JMxf/781B4AqhX4icZj3KpVq9C+fXuceeaZScfPOusstG/fHqtXr04YGzlyJDp16oT77rsPtm073sfYsWPx5z//GVdffTVOO+00rF+/HhdccIHnY9y6dSsuvfRSjB8/HmPGjMHixYsxduxY9O7dGyeeeCIA4L///S9ee+01jBw5Ejk5OdizZw+eeOIJDBgwAF9++WWFfRrp559/xoUXXohRo0Zh5MiRWLRoEUaNGoUXX3wRU6dOxcSJE3HFFVfgwQcfxKWXXooffvgBDRo0AAB8/PHH+Ne//oVRo0ahdevW2L59OxYtWoSBAwfiyy+/REZGBoCSF4VBgwbBMAzMmDED9erVw9NPP500Uv78889jzJgxGDJkCO6//34cOXIEixYtQv/+/bFp06ajnjQgImc5OTkYPXo0nnrqKUyfPt21jtx7772444478Otf/xrXXHMNfvrpJyxYsABnnXUWNm3ahIYNGwIoOYF25MgRXH/99WjcuDE++ugjLFiwAP/73/+wfPlyZZ+RSARDhgxB//798dBDD8VqxFtvvYWhQ4eid+/euPPOO2GaJpYsWYKzzz4b//znP9GnTx8AwMSJE7FixQpMmTIF3bp1w/79+/Hee+/hq6++Qq9evfD73/8eubm5+N///odHHnkEAFC/fv2jPi7vvfceXn31VUyaNAkNGjTAH//4R4wYMQLff/89GjduDADYs2cPTjvttNiJyaZNm2LNmjUYP3488vLyXKPaq1evxmWXXYaTTjoJc+bMwc8//4zx48ejVatWSecvW7YMhw4dwnXXXQfDMPDAAw/gkksuwX//+19GsolqodzcXOzbtw+2bWPv3r1YsGABDh8+jKuuuqpC9j937lyYpolbbrkFubm5eOCBB3DllVfiww8/jM3xWme9vrcrNWnSJDRt2hQzZ85Efn5+hXw9RFS3pVITf/zxR/Tp0wcHDx7Etddeiy5dumDnzp1YsWIFjhw5ktInvr/99ltcdtllmDhxIsaMGYMlS5Zg5MiRWLt2Lc4991zH261btw4XXnghWrRogZtuugnZ2dn46quv8Le//Q033XRTmY+DagmbjlkHDx60AdgXXXSR67xf/epXNgA7Ly/Ptm3bvvPOO20A9uWXX54wt3Ss1MaNG20A9tSpU5V5Y8eOtQHYd955Z2zbkiVLbAD2tm3bYtvatWtnA7Dffffd2La9e/fawWDQvvnmm2PbCgsL7Wg0qtzHtm3b7GAwaM+ePVvZBsBesmSJ69f89ttv2wDst99+O7ZtwIABNgB72bJlsW1ff/21DcA2TdP+4IMPYtv/8Y9/JNzPkSNHEu5nw4YNNgD7ueeei2274YYbbMMw7E2bNsW27d+/327UqJHy+Bw6dMhu2LChPWHCBGWfu3fvtrOyshK2E5E3pbXo448/tr/77jvb7/fbN954Y2x8wIAB9oknnhhb3759u+3z+ex7771X2c+WLVtsv9+vbE9WB+bMmWMbhmHv2LEjtm3MmDE2AHv69OnKXMuy7E6dOtlDhgyxLctS9puTk2Ofe+65sW1ZWVn25MmTXb/WCy64wG7Xrp3rHAmAnZaWZm/dujW27d///rcNwF6wYEFs2/jx4+0WLVrY+/btU24/atQoOysrK/Y4JKvJJ510kt26dWv70KFDsW3vvPOODUA51tLbNm7c2D5w4EBs+1/+8hcbgL1q1SrPXxcRVb/S2qv/CwaD9tKlS5W5Xt5DDhgwwB4wYEBsvfS9XdeuXe2ioqLY9vnz59sA7C1btti2XbY66/W9Xenx9e/f345EImV+bIjo2FOemjh69GjbNE37448/TthvaV3Tf2/X7zfZ7+SvvPJKbFtubq7dokUL+5RTTolt03+HjkQidk5Ojt2uXTv7559/TnocVDcxOn0MO3ToEADEPnHnpHQ8Ly9P2T5x4sSj3kdpnG7SpEnK9htuuMHzcXbr1k35xGXTpk3RuXNn/Pe//41tCwaDMM2Sp3M0GsX+/ftRv359dO7cORYTrAj169fHqFGjYuudO3dGw4YN0bVrV/Tt2ze2vXRZHmN6enpsORwOY//+/ejYsSMaNmyoHOPatWvRr18/9OzZM7atUaNGSmwSKPnr0MGDB3H55Zdj3759sX8+nw99+/bF22+/XWFfN9Gx6vjjj8fVV1+NJ598Ert27Uo659VXX4VlWfj1r3+t/CxmZ2ejU6dOys+irAP5+fnYt28fTj/9dNi2jU2bNiXs+/rrr1fWN2/ejG+//RZXXHEF9u/fH7uv/Px8nHPOOXj33XdhWRaAkmjNhx9+iB9//LEiHoqYwYMHo0OHDrH1k08+GZmZmbF6Z9s2XnnlFQwbNgy2bSuPyZAhQ5Cbm+tYl3/88Uds2bIFo0ePVj5dOWDAAJx00klJb3PZZZfhuOOOi62Xvl7I+ktEtcfChQuxbt06rFu3Di+88AIGDRqEa665JiGil6px48Ypn+bRa0ZZ6qzX93alJkyYUGP6kRNR7VDWmmhZFl577TUMGzZM6e1YKllc2ouWLVsqF4jNzMzE6NGjsWnTJuzevTvpbTZt2oRt27Zh6tSpsXRPeY+DagdGp49hpScQS084OnE6IZmTk3PU+9ixYwdM00yY27FjR8/H2bZt24Rtxx13HH7++efYumVZmD9/Ph577DFs27YN0Wg0NlYa5asIrVu3TiiKWVlZaNOmTcI2AMoxFhQUYM6cOViyZAl27typxM1zc3Njyzt27EC/fv0S7lt/zL799lsAwNlnn530WDMzM718SUR0FLfffjuef/55zJ07N2k/mW+//Ra2baNTp05Jby/ju99//z1mzpyJv/71r0p9ANQ6AAB+vz+hh1fpz/2YMWMcjzc3NxfHHXccHnjgAYwZMwZt2rRB7969cf7552P06NGuF/by4mg1+aeffsLBgwfx5JNP4sknn0y6j7179ybdvmPHDgDJXyM6duyY9Bd3/XhKTzrqjy8R1Q59+vRRfjm+/PLLccopp2DKlCm48MILy32Rl6PVjLLUWa/v7Up5ee9MRCSVtSb+9NNPyMvLQ/fu3Sv0ODp27Jjwe/AJJ5wAoKTndnZ2dsJtvvvuOwCo8GOhmo8nGo9hWVlZaNGiBT777DPXeZ999hlatWqVcOJK/hW3Mjn95Ve+mbvvvvtwxx134De/+Q3uvvtuNGrUCKZpYurUqbG/OlfmsXg5xhtuuAFLlizB1KlT0a9fP2RlZcEwDIwaNSqlYyy9zfPPP5+0sKd69W8iUh1//PG46qqr8OSTT2L69OkJ45ZlwTAMrFmzJmktKP1kXjQaxbnnnosDBw7g//2//4cuXbqgXr162LlzJ8aOHZtQB+QnteV9ASUX8pKfek52f7/+9a9x5plnYuXKlXj99dfx4IMP4v7778err76KoUOHlv2B+D9Hq3elx3jVVVc5/qJ+8sknp3z/ZT0eIqrdTNPEoEGDMH/+fHz77bex/typ8lrDvNTZsr63q6r3zkRUd1VUTXT6RKH8wA5Rqngm4hh34YUX4qmnnsJ7770Xu3K09M9//hPbt2/Hddddl9L+27VrB8uysG3bNuXTPsmu7lceK1aswKBBg/DMM88o2w8ePIgmTZpU6H2lasWKFRgzZgwefvjh2LbCwkIcPHhQmdeuXbukj4++rTS62KxZMwwePLgSjpiISt1+++144YUXcP/99yeMdejQAbZtIycnJ/aX3WS2bNmC//znP3j22WcxevTo2PZ169Z5Po7Sn/vMzExPP/ctWrTApEmTMGnSJOzduxe9evXCvffeGzvRWBmxlaZNm6JBgwaIRqNlrk3t2rUDkPw1oqJfN4io9ohEIgCAw4cPV/p9laXOen1vR0RUkdxqYtOmTZGZmYnPP//cdR+ln+Y+ePCgEmsuTZfotm7dCtu2lfeO//nPfwDA8QKkpfX0888/5++rxxj2aDzG3XrrrUhPT8d1112H/fv3K2MHDhzAxIkTkZGRgVtvvTWl/Q8ZMgQA8NhjjynbFyxYkNoBO/D5fAmfXlm+fDl27txZofdTHsmOccGCBQl/NRoyZAg2bNiAzZs3x7YdOHAAL774YsK8zMxM3HfffQiHwwn399NPP1Xg0RMd2zp06ICrrroKTzzxREIfmksuuQQ+nw+zZs1K+Bm3bTtWW0s/RSPn2LadNI7tpHfv3ujQoQMeeuihpG8uS3/uo9FoQmyvWbNmaNmyJYqKimLb6tWrlzTeVx4+nw8jRozAK6+8kvRNrlttatmyJbp3747nnntO+frWr1+PLVu2VOhxElHtEA6H8frrryMtLQ1du3at9PvzWmcB7+/tiIgqytFqommaGD58OFatWoVPPvkkYby0ZpWeBHz33XdjY/n5+Xj22WeT3u+PP/6IlStXxtbz8vLw3HPPoWfPnknTdQDQq1cv5OTkYN68eQl/gGHypG7jJxqPcZ06dcKzzz6LK6+8EieddBLGjx+PnJwcbN++Hc888wz27duHl156SWn8Xxa9e/fGiBEjMG/ePOzfvx+nnXYa1q9fH/vrR0V9mubCCy/E7NmzMW7cOJx++unYsmULXnzxxXL3IqtIF154IZ5//nlkZWWhW7du2LBhA954442EHpK/+93v8MILL+Dcc8/FDTfcgHr16uHpp59G27ZtceDAgdhjlpmZiUWLFuHqq69Gr169MGrUKDRt2hTff/89Vq9ejTPOOAOPPvpodXypRHXS73//ezz//PP45ptvlJhKhw4dcM8992DGjBnYvn07hg8fjgYNGmDbtm1YuXIlrr32Wtxyyy3o0qULOnTogFtuuQU7d+5EZmYmXnnllTL1EjRNE08//TSGDh2KE088EePGjUOrVq2wc+dOvP3228jMzMSqVatw6NAhtG7dGpdeeil69OiB+vXr44033sDHH3+sfPKmd+/eePnllzFt2jSceuqpqF+/PoYNG1bux2ru3Ll4++230bdvX0yYMAHdunXDgQMH8Omnn+KNN97AgQMHHG9733334aKLLsIZZ5yBcePG4eeff8ajjz6K7t27V8mnmYioeq1ZswZff/01gJJ+rsuWLcO3336L6dOnV0n/aa91FvD+3o6IKFWp1MT77rsPr7/+OgYMGIBrr70WXbt2xa5du7B8+XK89957aNiwIc477zy0bdsW48ePx6233gqfz4fFixfHfp/UnXDCCRg/fjw+/vhjNG/eHIsXL8aePXuwZMkSx2M3TROLFi3CsGHD0LNnT4wbNw4tWrTA119/jS+++AL/+Mc/KuZBohqHJxoJI0eORJcuXTBnzpzYycXGjRtj0KBBuO2228rdvPW5555DdnY2XnrpJaxcuRKDBw/Gyy+/jM6dOyMUClXI13DbbbchPz8fy5Ytw8svv4xevXph9erVSfupVZf58+fD5/PhxRdfRGFhIc444wy88cYbsU99lmrTpg3efvtt3HjjjbjvvvvQtGlTTJ48GfXq1cONN96oPGZXXHEFWrZsiblz5+LBBx9EUVERWrVqhTPPPBPjxo2r6i+RqE7r2LEjrrrqqqR/6Z0+fTpOOOEEPPLII5g1axaAkp/l8847D7/61a8AlFwUZtWqVbjxxhsxZ84chEIhXHzxxZgyZQp69Ojh+TgGDhyIDRs24O6778ajjz6Kw4cPIzs7G3379o21ucjIyMCkSZPw+uuvx66K3bFjRzz22GPKlawnTZqEzZs3Y8mSJXjkkUfQrl27CjnR2Lx5c3z00UeYPXs2Xn31VTz22GNo3LgxTjzxxKTxc2nYsGF46aWXcNddd2H69Ono1KkTli5dimeffRZffPFFuY+NiGq2mTNnxpZDoRC6dOmCRYsWpdzGJxVe6izg/b0dEVGqUqmJrVq1wocffog77rgDL774IvLy8tCqVSsMHToUGRkZAErel65cuRKTJk3CHXfcgezsbEydOhXHHXdc0t8jO3XqhAULFuDWW2/FN998g5ycHLz88stHrXdDhgzB22+/jVmzZuHhhx+GZVno0KEDJkyYkOIjQrWBYfMzq1QNNm/ejFNOOQUvvPACrrzyyuo+nFph6tSpeOKJJ3D48GHHRuZERHVVz5490bRp0zL1tCQiIiIioqrFHo1U6QoKChK2zZs3D6Zp4qyzzqqGI6r59Mds//79eP7559G/f3+eZCSiOi0cDseanJd655138O9//xsDBw6snoMiIiIiIiJPGJ2mSvfAAw9g48aNGDRoEPx+P9asWYM1a9bg2muvRZs2bar78Gqkfv36YeDAgejatSv27NmDZ555Bnl5ebjjjjuq+9CIiCrVzp07MXjwYFx11VVo2bIlvv76azz++OPIzs7GxIkTq/vwiIiIiIjIBaPTVOnWrVuHWbNm4csvv8Thw4fRtm1bXH311fj9738Pv5/nupO57bbbsGLFCvzvf/+DYRjo1asX7rzzTgwePLi6D42IqFLl5ubi2muvxfvvv4+ffvoJ9erVwznnnIO5c+emfGEyIiIiIiKqGpV2ovHAgQO44YYbsGrVKpimiREjRmD+/PmoX7++420GDhyI9evXK9uuu+46PP7445VxiERERERERERERFRBKu1E49ChQ7Fr1y488cQTCIfDGDduHE499VQsW7bM8TYDBw7ECSecgNmzZ8e2ZWRkOF62nYiIiIiIiIiIiGqGSrkYzFdffYW1a9fi6aefRt++fdG/f38sWLAAf/rTn/Djjz+63jYjIwPZ2dmxfzzJSETHsoULF6J9+/YIhULo27cvPvroo+o+JCKiasF6SETEWkhENV+lfKJx8eLFuPnmm/Hzzz/HtkUiEYRCISxfvhwXX3xx0tsNHDgQX3zxBWzbRnZ2NoYNG4Y77rgDGRkZjvdVVFSEoqKi2LplWThw4AAaN24MwzAq7osiojrLtm0cOnQILVu2hGlWyt9fUvLyyy9j9OjRePzxx9G3b1/MmzcPy5cvxzfffINmzZod9faWZeHHH39EgwYNWA+J6Khqai0EylcPWQuJqKxqaj3ke0Miqkop10K7Etx77732CSeckLC9adOm9mOPPeZ4uyeeeMJeu3at/dlnn9kvvPCC3apVK/viiy92va8777zTBsB//Md//Ffufz/88EO5619F6tOnjz158uTYejQatVu2bGnPmTMn6fzCwkI7Nzc39u/LL7+s9seU//iP/2rfv5pWC227bPWQtZD/+I//KupfTauHfG/If/zHf9Xxr6y1sEyX/J0+fTruv/9+1zlfffVVWXapuPbaa2PLJ510Elq0aIFzzjkH3333neOVJmfMmIFp06bF1nNzc9G2bVv0uOR2+AKhlI+lrmj0/s7qPoQa5cAZrar7EGqMN+9eXN2HUGPkHbbQrtd2NGjQoLoPJaa4uBgbN27EjBkzYttM08TgwYOxYcOGpLeZM2cOZs2albCd9bAE66GK9TCO9bBETayFQNnr4dFqoWFpA2LdjNqxZSOqTjOs+JgZsdWxqFyO71C/L1OOReR9qRPlfSHhvqJiWc7TDliZJ8a0+4Iv/gkF2+fTxsS63yfmqZ+EUm7nj4/Zpj5P3JdfHbPkmPjQhLxNybq4jet9yX1rn9yS+xfLkXR1XiRkiLH49mhQ3Z2VHv8+REPq98sKiu95eiS2HAhGlHnpoXBsuX4wntbKChUq8xql5ceWG4tlAGiadii23MyfG1/25SnzmvsPx5d98edGU189UImaWA8r8r1h98vugC8tVHIKQZA1Sy6bljZRlBRZNwG1HspaaWi1TO5Tjpl6PXSrlXJd1mFLq3MyxKm/BnglS5H8NKj2yVBbrsvb6DXKdKlfcswva7RWN2UNFGNWwHmepddeORaIb4+mQZ0n9hkVv1Yk1MOgcz2MhsSDnx7/hgVE/QOAUCheHzNFDdTrYZO0eC1rHFTrYeNAfKyJP14bjzPVeVm+gviyGd9/PUOr0eJhCxnq65JfrAeM+IuPmWKnQks8SaNaAFkZEz/AtjYvLMbkPrR3CsqPg/ajDAA4fNjCmX33lbkWlulE480334yxY8e6zjn++OORnZ2NvXv3KtsjkQgOHDiA7Oxsz/fXt29fAMDWrVsdTzQGg0EEg8GE7b5AqKR4HuP8ZuJjcyzjcyIus0HNiYHUFDUpQrJv3z5Eo1E0b95c2d68eXN8/fXXSW+j/+ElLy8Pbdq0YT38P6yHKj4n4lgPVTWpFgJlr4dOtdBukA47LZTwLtuU6y4nGpVfmF3Hki8DgBkWv0z7k590LLmdGDP1X7rFmPxibP3XB3GiUY5Z+hlU8fw31BONtinWTXmiUT/5J+cl/wVZv53riUavJxAdfsku2Z84JI8nIfWXCVOcaPTJMb18il+mbe0Xa4hfrI2Q+J5ov1gj5BNj8fs1tGMyxIlLpGnPjbT4L8ZGIL5/n79YmWf6wmI5/ku236f+Yn0s8//fiaqaVA8r8r2hUS8EI839Dy/KSUe9bIjaZmpPG1MZkycTtXny5KIha55WD8WYqR2wYci5Yh/6FyZPwshzhC5d5Gy3773TSUdtXalLpl4PK+BEY0COye3OJxOj2pgh18WJRmUZAIJinqxL+ttrl3pop4vvi6iNZrp6WspMj9csU9Rbn7Y/Q/zBxpemPsF8ATEm6mFAq3NpYt0v3hCkad/WgHhOBRLG7KTzfPqZfEGehLRczn5HtX3Ic/7yVUQ/SWiKJ7ollsPavKgY07/lUdtAwE6tFpbpRGPTpk3RtGnTo87r168fDh48iI0bN6J3794AgLfeeguWZcVOHnqxefNmAECLFi3KcphERMckpz+8EBEdS1gLiYhKsB4SUXWolD/hd+3aFb/85S8xYcIEfPTRR3j//fcxZcoUjBo1Ci1btgQA7Ny5E126dIldJeu7777D3XffjY0bN2L79u3461//itGjR+Oss87CySefXBmHSURUYzVp0gQ+nw979uxRtu/Zs6dMnwwnIqrtWA+JiFgLiaj2qLSs0IsvvoguXbrgnHPOwfnnn4/+/fvjySefjI2Hw2HqViuyAAAgAElEQVR88803OHLkCAAgLS0Nb7zxBs477zx06dIFN998M0aMGIFVq1ZV1iESEdVYaWlp6N27N958883YNsuy8Oabb6Jfv37VeGRERFWL9ZCIiLWQiGqPMkWny6JRo0ZYtmyZ43j79u2VhpVt2rTB+vXrK+twiIhqnWnTpmHMmDH4xS9+gT59+mDevHnIz8/HuHHjqvvQiIiqVEXUw0jIgB00Ei7kYot2TUbUSLodUC8aovcrk2OyX5nWdkvrwmSKJb2fmOjR53ahAbFs6Hcmj1H2IXPpSZYw5nQ77XBlTzUbFdvTzjb09RT2n7APh7FUD912WNbWlYfT0vq1iYOKirGwpfbNLBbrRZb6q1yhuJLDEStNLKvR2XwzfrGZI3a8F9phW+sbeQzLt1O9Ykjlqqj3hlbAKOnNp18MRvaYlQ9BQhtCuUF/0htJh0ytvsiHWPYeNPQLz8ieh9rPv+HSA08h+0261UC5b7f+jZbsveiyPznP0ObJx1ev3ymwXR4nvY56UhGlXP9WigNRyqZ2gEo9FMsRrW5GRD0M22qtlOtFojYWmmonwkI7vh4Q/YzTXK8apI7JPopu/RbVo9f7Kiffu34xmKjjPOd9FNuyH6T6WMsejVaSJ0pRSk+eSjzRSERE5XPZZZfhp59+wsyZM7F792707NkTa9euTWgCTkRU17EeEhGxFhJR7cATjURENdiUKVMwZcqU6j4MIqJqx3pIRMRaSEQ1H080EhEREVGdFw0BCAJ2RI0BWTLqHBYRYK2TuS2ibXr8GkqcGcmXocamlFSiHk0S9yUjhYAWZ5Q7dE5hqSwt1iWjaJXWvb2SyISmnu5SYoRluJ2H+3LlEp2WMUp9miXGZHRajwoWR+O/vhVr0WkZpZZRwXwtOp0h1g8Z8eh0yCgAlTisx3frmGgagKBWTwDlOWqKmmJr8WDZIgJ6OwMlIuz8gyPvWz7csoUFoLZmcGuDYDgVWJ3X9hEubRqU4Lil128RPxf7S7zb5BFzNymmWJ1fo/R9VmznC3d20sWSdYdWElHthTki1vVWEnK90BbLlhadNuJtJgLihdRn6vFoETHWXnAD4glXJL7R+kuqz+PjK2PQ+lNZHoeMlbtFosPicYpq8yz9zY62jyMptpGobW8niIiIiIiIiIiIqAbiiUYiIiIiIiIiIiIqN0aniYiIiKjOiwZREp1WL0wJQ6zbSmRZnSejg7Z2hVDPcTblCqzxZUvLU5kiymRrEU7D4arTrsQ+bC06bbhdkVq50rSIA2oxStvrVa0rUxnigLZDjNAtYu1Kpi3drrKqXHVai6/JeKBDjBrQooJR9Ve5gqi40rTP+arTcl1egTqDV52OqalXna4opfVQz6wqcWmxbGg1SqYt9Vi17XhFai0S7VQP/Xq7CHFF6lSuOg+XK0jr25V6qz84ye9b32o7XWla+5mXUe8Kr5ouD5PX1yv9ytUptY/Q78xO/ngk1kNZA52vOu3WSuKIFa+BwWg86x8yIso8X0L/gP87Bv0a0aJWhm0tOm1ExXJ8fz79O+vxG320K0GXKhbHqEeglYi113lJYtT5lte+LCp+opGIiIiIiIiIiIjKjScaiYiIiIiIiIiIqNx4opGIiIiIiIiIiIjKjT0aiYiIiKjOs4I2ELJhh7V+iPLP7k7LgNoPUW/Ep050HpK9/ERbKNuv3pnso5jQD9J06Fdm6H3NHPomWi59GLV+ZXIfeh82T/T78iWfViZKT0VvTcPcei+69XJ03qHHefpc22tPsvjzIaz1JAtH4w9isU99QItEj7JCKxBblr3KACBf9GgMWenxZYM9GkvlW3W8R2PIhh2yYUa1+iL7MopWdqb+s6vUyoQuhWJZ/rDp9UUsi7MStn5Msneudl+GXPfaSk4eR7QM3+eoLODifhOOSfTYVfq3Ove2TeghKWuFyyF5rV+ufRlTqIGe+xIn3DD5sls9lP1rZf0DgGIrvq73rC32xddlDQxZap0zxRPRhPPzQfYyDGhPtjSlR6PaA1JK6Nno4b6i2psR2WPRrb+i7DHpeV6SHo1H2KORiIiIiIiIiIiIqgtPNBIREREREREREVG5MTpNRERERHVeNAjYQcD2afE1JZYntmvRMHXdJSqoROXUWYaMg4l34UZUPyZ5I+2+EmKKZWRr0TAZEdViunCKSyfE/DxmiW2HSKWbVKN8XmOEXiOFXukPhUNUUJ9nO0Sno9r3REYFiy31V7kiMVYQFdFpU41OH1Gi0/EYYchIB5VINS5YW0RDgB0CLK29gUiAwhQJUL3lhFIrEz66JOfKeLDe3sFhWYtp2z4ZU9bGRH1U9u5Wo5zaSuj0MadWDXraVhZ+GQnWXxC81s0UeI5Ka3MrupVEQpcRj/VQaSUhDiqit5IQNa9Qi07LGugXT+xDRsj5gOW+bXV/ITu+P73NhIxL+8ST2WtUWhd1efBlvNkt9iz3ocattXni8bWSfA6xwHaOgrvhJxqJiIiIiIiIiIio3HiikYiIiIiIiIiIiMqNJxqJiIiIiIiIiIio3NijkYiIiIjqPCtkAekW7GKt15jsT+XWG9G1X5XciVtPMjEm+nopPci09YSxiBhz6k+ms116krmNpTJP9ivzOc5y5dpfTEq1D6NDTzLbqQfb0Q5D6TWmN0BL3q9N9mQEAEv2yRJjkajWkywqejRG1Qe4UPQkKzDjPcmCptpjS/YoC5qyR6Pad+xYVtd7NFrpJfUQ2pdpyPoSlnXIpbetx9po2Ho/yPiYeLoepY+uVlMNpxUXLrXMsNx66sXHlONwuV9DjLm2ZNTvN4XamdgrUx6ImFem17YU5rn0XnRqCJlQD2WfWlEDw1o9LIrET2cFzIAylhaNP6n84gXXTGgcGSf7HIZ9RcpY2Bb9IA31vtQejfH9+xIaeJad156KiT0aRS9HcRt9f5bLEycKA4VR9mgkIiIiIiIiIiKiasITjURERERERERERFRujE4TERERUZ1nBy3YQQswtNiQY/ROj1cZjkNKEkvsUI/hGXJMRAUN7R25El80XWJublFvJ9ox2UrUO8VYdQXEpb1yTHm5fP16VNAxVq2nnlNIUiek8myH55S2cxkdlNHpqOUcmwtb6oNdHI0/kYp88bEiS32CHbHSYsuHo/EYdZqRWkSuLiqw6/ZjYWdEYKdHgIj6/LJFnFm2mUhoOeG19jjURgAwovFBmVQ3I/p9ichyQuy3fHHphKi0x5on22LYPq0oKe0SRHE0tOIohxJ24ZazlgfibbtrKwmn26XWScKV4fTwWs71UEanI1q7iLB44hRFtToXidc5M+H1PDkZMZZRaQAoEm0mAobac0CumyKmXRnRaRl1jopvkh6BdhrTI9aWyzc6apsoiqTWUoOfaCQiIiIiIiIiIqJy44lGIiIiIiIiIiIiKjdGp4mIiIiozjODEZihCCwtvmYrmTUlD6jtwTn2qkanxXYtNSXTVqY4jIQ4sOk85vVqpwoZ37O1g7IcYn6Ac3TQ4xVS9RixxzCgthPnIddoc1VGBV2+MOeooPM+3KLTMjoYNtWdFIu4tNMVqAE15nfEjMcLA9H0hOM/VqV6pdXaIhAKw0z3wdKiqFY4/nyzTXHVWi0e7HqFY/GDKWugljaFJR5iyy9intpVp03lqtOOd+XOcvhB9FrzXBhRdR/KHkSNNkznmLbbvaotLVy+4OQXdD7qmMOFoF3nuVHqfkIrCTFP1Dn9qtNKKwkRndbrYXEk/vz1aS84ypWmPb76RMQLc7qvWBmTLSgCLjVV3pdPfxPgkR5vlmTUWc7TI9DKmMtVp5V9J/kmF0VT+2xilXyiceHChWjfvj1CoRD69u2Ljz76yHX+8uXL0aVLF4RCIZx00kn4+9//XhWHSURERERERERERCmq9BONL7/8MqZNm4Y777wTn376KXr06IEhQ4Zg7969Sef/61//wuWXX47x48dj06ZNGD58OIYPH47PP/+8sg+ViKjKvPvuuxg2bBhatmwJwzDw2muvKeO2bWPmzJlo0aIF0tPTMXjwYHz77bfVdLRERJWDtZCIqATrIRHVFZV+ovEPf/gDJkyYgHHjxqFbt254/PHHkZGRgcWLFyedP3/+fPzyl7/Erbfeiq5du+Luu+9Gr1698OijjyadX1RUhLy8POUfEVFNl5+fjx49emDhwoVJxx944AH88Y9/xOOPP44PP/wQ9erVw5AhQ1BYWFjFR0pEVHlYC4mISrAeElFdUak9GouLi7Fx40bMmDEjts00TQwePBgbNmxIepsNGzZg2rRpyrYhQ4Yk/EWn1Jw5czBr1qyKO2gioiowdOhQDB06NOmYbduYN28ebr/9dlx00UUAgOeeew7NmzfHa6+9hlGjRlXloRIRVZqqrIX+YAS+YARhrQWRbA2ldlNybgaW0HvRoQ+ZHdX6TkXjvZusSHzM1Hocyn5lpn7Asi+j1x6NSk8yvU+Y/GL0A4nPNcSy7dCTEVD7iaXUk/FoHL5kr33H9H147TvmuRec155ktnNPMlv0IbO0HnIR0S8rrPXNC1vxb0xxNP5rXoEZUObJ/mL+aNn7mB0LCqPhKr/PqqyHGelh+NJNhLUejWF/fD0ilq1C5962htZPTq2Hom7qPRpFbZN9GW2f9rMhezaaLj+IhmMxV2qg4dSvEQCiHnvqyfvy6V+/+FpMl/uyXcbKS3ttUMqNSz303KfWa79GfV22C3bpWStfO+XrYVTrFyh7KhZrPTD9ZrwGGnrDYOWu4/uPiO9lWHteB8TB+916NLrcly/hiZmc1z6Kbr0X1Xleez4mfmOLIzWwR+O+ffsQjUbRvHlzZXvz5s2xe/fupLfZvXt3mebPmDEDubm5sX8//PBDxRw8EVE12bZtG3bv3o3BgwfHtmVlZaFv376Of6QB+AlvIqpbWAuJiEqwHhJRbVIlF4OpTMFgEJmZmco/IqLarPQPK2X5owtQ8gnvrKys2L82bdpU6nESEVUm1kIiohKsh0RUm1RqdLpJkybw+XzYs2ePsn3Pnj3Izs5Oepvs7OwyzSciohIzZsxQWk/k5eXxDSURHXOcamFGKAxfyESBlgwqFstqok6PtiafB6iRWKcYNQAYkfiyjEtbWhTZFPFAW48KOnxMwNaick7JNluL6Ble43su8yo9Lu0klcgfnGOEbhFr9x06LLuN6VFBh1i1ZakHIaNtES32WhSJ/2rnc4v5iSi1XzxJA/oT9hhWZJXhCVWDOdXDzFAh/Ok2irTnUKEv/twoFM8TvYODfKboz1HZj8IQ7SIMrZWEJRL9tkiq6/XQdquHXttHKDt0ayUhx1xirqbL57XkMbnVV9d663S/+tcvbuL2WHhuJWEk3Z50PRVONVA/KFkPxfMmqn39EfF9MA3tuewQYdZf25XotGhbkWaqp8pkHQ1o/VPkmM8lOu3UnsIqw4MbdYhEWy5fl9Pt9X0kHJdtIFwTo9NpaWno3bs33nzzzdg2y7Lw5ptvol+/fklv069fP2U+AKxbt85xPhFRXVP6h5Wy/tGFn/AmorqEtZCIqATrIRHVJpUenZ42bRqeeuopPPvss/jqq69w/fXXIz8/H+PGjQMAjB49WrlYzE033YS1a9fi4Ycfxtdff4277roLn3zyCaZMmVLZh0pEVCPk5OQgOztb+aNLXl4ePvzwQ/7RhYiOGayFREQlWA+JqDap1Og0AFx22WX46aefMHPmTOzevRs9e/bE2rVrY/0lvv/+e5jiI6+nn346li1bhttvvx233XYbOnXqhNdeew3du3ev7EMlIqoyhw8fxtatW2Pr27Ztw+bNm9GoUSO0bdsWU6dOxT333INOnTohJycHd9xxB1q2bInhw4dX41ETEVWsqqyF9dKK4Q8mRolkei2sJOq0eJGIVEGLChoBEQ8UV5M2tHfaprh6puWP38bULnArk0x6qklG22yXlJtyG7c4oLxCqlusurbxeJVV16uxOu26DA+L4ZScdIsKypi+dpVVedXViHYg8irU8grUhVH1qtN+M/mVpt2ulnqsKY5W/aUMqrIeHhc6gkAogiORNGV7QD43xHKB9twoUkqKVivlFYNFPUyoc6I+yri021WnE2K/Xr9NThHmhHrodFlkOM/TI8vKPuTrhsfWFKlyiUe7RaI9t5JIhduXJe9A//qV6LSIB/vUeVERzY9oV502IvEnleerTovzUhGt5YSMR8saCgBm1FsdNQ1vV50+Wpw5tuzyouV0dWr9vc3R7jccSTLRg0o/0QgAU6ZMcfxE4jvvvJOwbeTIkRg5cmQlHxURUfX55JNPMGjQoNh6af+cMWPGYOnSpfjd736H/Px8XHvttTh48CD69++PtWvXIhQKVdchExFVONZCIqISrIdEVFdUyYlGIiJSDRw4MOGTI5JhGJg9ezZmz55dhUdFRFS1WAuJiEqwHhJRXVH1nwknIiIiIiIiIiKiOoefaCQiIiKiOq9BsAiBoJ3Qn8gW63I5rPUds6OyZ5T6qSPZk8wU7fBstcWT0ofMdOtJJntN6f20TLHBdGmi5fTJKH270hBQH3Poa5awD6/zxLLeqkr+VqL0odQeG4celW49yRJaUnnsy5hKjzL9JmpfRodloAw9yWQPMXUsEo0/qYoM555hfoM9Go+mOJLCN78WaZR2BGnBCEI+tQFbmi9etOTzwU7oWRtfL46on12StdIKi+eyX9uH+JmXNdDWntdufVT1+uCJS49G23Kph/IwRD9URLVCL47JkD1wfeq0hBrowO1rdKxzXvvS6nNdH2uHfZTlWyBuqPQA1l9v5eOm1EbtOWTInrUeD0F7ACKi/3JEfF+jpvq89okXbb1WynW1pnr8Jrtw7dcoHvyEvtIOL2Bu74GSiURS+xr4iUYiIiIiIiIiIiIqN55oJCIiIiIiIiIionJjdJqIiIiI6ryGaQUIpEUTY0NiOSrGrKj69/hINB4f0tJ2MGTUVSQRDS0SbYp33pYY0+OxMh5t6/FoJb4mY9RwZjlHBZUxWx8rf+yr2niOALpErL3elUtaXIkKyodTe2hTiQpGtYhpsRGP9hku0WmfjE4ry4xOlwq7xGbrgiZphxFMCyBopivb5XNA1koZL9XXIxE1ExwVcWlbxqP9LnVOtpLQaplSAz1GghN+lJWYrssPrFs9lLsTaWklRq3v062VRApcY88u89zj597mlSkiHduhtgunh8Olbsrvia29LitPS62liWF4+0ydfC1W3g9YenTauVY6RaeNFGuqW5zZcvhGeI1Ou+47yRij00RERERERERERFRteKKRiIiIiIiIiIiIyo0nGomIiIiIiIiIiKjc2KORiIiIiOq84wJHkJYWTuhBJPswhaPxRmHRoNYLSvTNsyJaL6SA6MkUjo+Z2jttS6ybbj3JfMmXgaP00HIieo3Zep8wsW649RBz6zVWAb3HUmK4NWWLc+1rVtN7klnOfbeiWr8yU/RsjIixsKnOK4wEkhw4SeFI3e7R2DBwBKGA37Uvp3yuydoIAMViPRxQC50ViI/ZftG/VutZ61TnEmqeKZe9/VDqtcxw6MuYWA+tpPMSyZqqHpPhsUejnFchzzavtUzv0ejQA9NrP8iUudZDsSx71upthGUPW60nYzQKT+T3Lyruy6f1wJU9Gn0uPRq9bC8Lp16LRxuzvc47SnPiqPaz7xU/0UhERERERERERETlxhONREREREREREREVG6MThMRERFRndck7RBCaQFYWk652PIlXY5osdRIREQF09QxW8S3bBGjtsLqMZgiopZyVNBjdNCRpUW5LMt5TEb99DEnlR2jdvjy9fSXWxospfi5G2WH+mPosKxFAJ2jglrcVDzXLO0zI1ER9Qsb8SeV4RLfcxs7lkXqeHQ601eAdJ8fPu35KiOWRaLXQ6EWjy4Q8fsCv5pRLfbHn9yWX8SD9Trn0D6iLJFd20jhB9itrnmtefJHz21eRdRDpV5p9cDhcXOrh26Pr1ttdBo7SvJW24nYhUs3Ds9106VW6vXR8ZDEF+ATy7aW03aLVTtFpCuivrpFm71Got2fooxOExERERERERERUQ3FE41ERERERERERERUboxOExEREVGd19h/GOl+Pwot9Yq7Mh5YLJe1uJByReqIOhYNiyuryitLa++0bbHuHp0W8S0teqVeIdRhWadc+dRyHku4mbgitdhuaDks2ymXVdGxQX33bokvj1dgregYtZ6Ukw+BoVzQVrtSrRyTQ1F1ni2urGoZWrRPzDXEPBmjdnO0CN2xJBLRM5p1SwOzEBmmLyF6Wd8Xr49H/Gmx5fxImjIv6I/ElgM+NTotr35uyYip6VzL1HYR6rEqsd/K/piU5VIrFS4/U7Idhds8N15/FJ1i1WW46rTTWKVfdVpKuHK3GHK5OrXyLdLbTMivRYxFEz5rJ3ci6qt2TDIebWmveTIiXdntKJzqdMoR66PcX9RK7YeOn2gkIiIiIiIiIiKicuOJRiIiIiIiIiIiIio3nmgkIiIiIiIiIiKicmOPRiIiIiKq85r4DiPD70Oh7dyjsSAaHyuKqG+Ti8V6OKD23bIC8b/d2wHRd8ynznPqy6j0XQRgi15mrv3KUuiZpfdTNGQ/Mb0nmeXQKKsiei+myPFr9tqv0WUs1Z5kjv3EtF3IHlqJfS7Fivye633H5Dytf6MlejFGQOURjabYW6+WyDCKkGH6UKw1iD1iBmPLQTP+LAr61GeUT/QH9Wm9Fw1TPEndapm8Xar9Bb2SP5huNc8r+fNbyR/dssVjo9copzHbcOlXmLAPseLyWDvuw+v+ysLpJUbfLtcTejSKnori8TC1b7ns2Sj7Epv6RNlHVDsMw+GFyUyxX6NbT0XJbe+p9HJM9tIejbJHIxEREREREREREVUTnmgkIiIiIiIiIiKicquS6PTChQvx4IMPYvfu3ejRowcWLFiAPn36JJ27dOlSjBs3TtkWDAZRWFhYFYdKRERERHVQE/8h1PObCdHpI9F4VLAgkBbfHlHnFYj1woD6FjoSiEesor743/GtgJpDEiltWD4ZsVaPVUYM9Vi1Gm3zmEtTItBqHExGxfQ4rzbR4315m+bGcyTcZZ57VNAhiuj1fsuShpPpUCXmp00TsT9DJjv1Y4rKXaufGbHkDY34WMQlRy1jdG6RumNNNBI9+qRaLM2IIs0o+V8KiHW5bGpPer+IlRpaPNTwGMWVauRTT6+Hei0uJyXerNdypyi5x1i5HlNXx/RWHQ6389qOoiwtJ1J5CD3W24SXKNmqQr4EutxQTdjrD6J8zqtfiP4z4HhfKXCPOpctBn2025RlP0dT6Z9ofPnllzFt2jTceeed+PTTT9GjRw8MGTIEe/fudbxNZmYmdu3aFfu3Y8eOyj5MIiIiIiIiIiIiKodKP9H4hz/8ARMmTMC4cePQrVs3PP7448jIyMDixYsdb2MYBrKzs2P/mjdvXtmHSURUZebMmYNTTz0VDRo0QLNmzTB8+HB88803ypzCwkJMnjwZjRs3Rv369TFixAjs2bOnmo6YiKhysB4SEbEWElHdUqnR6eLiYmzcuBEzZsyIbTNNE4MHD8aGDRscb3f48GG0a9cOlmWhV69euO+++3DiiScmnVtUVISioqLYel5eHgDgiZkLUL8BW1De0O6M6j6EGqWgSdvqPoQao/sHV1b3IdQY0SNFAOZW2f2tX78ekydPxqmnnopIJILbbrsN5513Hr788kvUq1cPAPDb3/4Wq1evxvLly5GVlYUpU6bgkksuwfvvv1/m+2M9LMF6qGI9jGM9LFHVtRCo2nrY3JeP+j4ThZYaic4PxKPTR6x4dDo/mqbMOxKJrxcE1H0U++NvqaPiitS2X7uysD95XDohOq1ckVobk7E3JebmMQ6l56AsGRvTIpBOV2TV9+ExW2XIiJqnWySJADrNS4geuh2It3kpxTndrorqtAwtLi3TrAnH4BKVE4+BJZ8cPvU2MknN6HRyVR2drur3hhVJv7KuU4y01qmIqLRbPNpFKld4dos9u0aiHe4roRw43JfbPM9XrtZVxJWr5f3Kp6T+UqHUQLFdex6r9THhutNJb1cRNdXrPtxehlONX5fl/nWVeqJx3759iEajCZ9IbN68Ob7++uukt+ncuTMWL16Mk08+Gbm5uXjooYdw+umn44svvkDr1q0T5s+ZMwezZs2qlOMnIqoMa9euVdaXLl2KZs2aYePGjTjrrLOQm5uLZ555BsuWLcPZZ58NAFiyZAm6du2KDz74AKeddlp1HDYRUYVjPSQiYi0korqlxn3EpV+/fhg9ejR69uyJAQMG4NVXX0XTpk3xxBNPJJ0/Y8YM5Obmxv798MMPVXzERETlk5ubCwBo1KgRAGDjxo0Ih8MYPHhwbE6XLl3Qtm1b10+DFxUVIS8vT/lHRFSbVEQ9ZC0kotqO7w2JqDar1BONTZo0gc/nS+gdsWfPHmRnZ3vaRyAQwCmnnIKtW7cmHQ8Gg8jMzFT+ERHVFpZlYerUqTjjjDPQvXt3AMDu3buRlpaGhg0bKnObN2+O3bt3O+5rzpw5yMrKiv1r06ZNpR47EVFFqqh6yFpIRLUZ3xsSUW1XqdHptLQ09O7dG2+++SaGDx8OoKRwvvnmm5gyZYqnfUSjUWzZsgXnn39+ZR4qEVG1mDx5Mj7//HO899575d7XjBkzMG3atNh6Xl4e31ASUa1RUfXQqRbuitRHvYgPP0aOU+bvDcf/SP1Tcf3Y8s9FGcq8Q8XxXo4FxWqPxkhYNFUsjv8d3wirvY2MsFgWLeCMCNR5ckxrk2iInoqQYx77JOp9wgyfmXS5ZINYl7cztXk+h36A2n3J/pK21jdQWffYa8xpOdUxz/P0npqm7TLmcR8+sQ/Z21P/WIhffNN9Wm88sW6KeaapPon8Yszniy8HfFXbl7Ami/ojR59USarivWGx7YPf9qFYeyKGxbpctrQfxIgVf2JGLfVJalmyER+SL5ztyLgAACAASURBVAOQLfAMl3mKmtL+0a1/o14fS+k9Gt167Cr1VtRNl160XvsruvVodOrDmLBPZZ7tOC/xOOQ3XWxP2IfDN9rlYU98CO2kY4ltb+U85yeY0/7Kwmn/bn0Q3Y5J3k4/Jqd+k/p9ed1/WVXqiUYAmDZtGsaMGYNf/OIX6NOnD+bNm4f8/HyMGzcOADB69Gi0atUKc+bMAQDMnj0bp512Gjp27IiDBw/iwQcfxI4dO3DNNddU9qESEVWpKVOm4G9/+xveffddpQdtdnY2iouLcfDgQeUv10f7NHgwGEQwGHQcJyKqqSqyHrIWElFtxfeGRFQXVHqPxssuuwwPPfQQZs6ciZ49e2Lz5s1Yu3Zt7AIx33//PXbt2hWb//PPP2PChAno2rUrzj//fOTl5eFf//oXunXrVtmHSkRUJWzbxpQpU7By5Uq89dZbyMnJUcZ79+6NQCCAN998M7btm2++wffff49+/fpV9eESEVUa1kMiItZCIqpbKv0TjUDJX2acotLvvPOOsv7II4/gkUceqYKjIiKqHpMnT8ayZcvwl7/8BQ0aNIj11snKykJ6ejqysrIwfvx4TJs2DY0aNUJmZiZuuOEG9OvXj1cVJKI6pSrr4a5oQ6RH/NgTzlK27y1uEFv+uTgel84tDinzjhSlxZaLi9S30FZxPGJoROJRI1OLTptyTCQzDS2xqkantViTEjH0mCOUMT9D+5yBjPlpY0qUWizbemxQyaUZjvOU6LTbmOd9IOly4pi+D4dlLc6sxAPFWMJ9ydu5jCnxaC32DKd9+LXsvEM8GlCj0zISLZcBwO+PP8H8Ilad5md0ulS0imPkVf3e8IgdBCwfjljqpx2PWPE6VxANJF0GgLAVf8JGonp0WqxHRbQzoQ2E+AETY3qSU94uYR+iBhqp1May0GtnKY/10G2e3kpC1gDboQ7pY27xaGUfZWgz4XhMLu0tUml94R7ndsrYQ41cm/oTRyzK2LPWSsIUt5PLeqTYcNhfwv5Rdm7xZZ2c6R6JRtIx96h04nHp0XivquREIxERxS1atAgAMHDgQGX7kiVLMHbsWAAlf3QxTRMjRoxAUVERhgwZgscee6yKj5SIqHKxHhIRsRYSUd3CE41ERFXM9vBX1lAohIULF2LhwoVVcERERNWD9ZCIiLWQiOqWSu/RSERERERERERERHUfP9FIRERERHXernBDhMJ+7CpWezTuK6ofW84tSo8t5xenKfMKi+M9yiJhrVFWcfxv90Zx8j6MAGCIdTlmRrV+T5Yc0/ZhOfQkcyOaNRk+7djluudeY+pnFWyxLnuN2T7tMw2yD5lbDzG3fogOvRfd5iV8tMJhzHbpJ6bM0/oruvaKVPoyigG915rcp9My1L6Mhjbm1JfRr/VeDIj+gwExL62K+xLWZJE6/ljkRtNRHPUjN5qhbD8UjfemzY/E+zceiWj1MBI/jVAcUZ/MVsShR2NE6yEnHmJZ5xJ61sqWenrNK29fxoS+i5bLmLyd7HnnsW66zLO1MdtxH+ounHvbus3Txhx6O7rXXtt5npF8HgCt9rr0rDW91UNZA019zJBj8e+rqR2TaSYf03sZ+gznMaf71ekvsaX0Vsxu9F6MZZ1nudw+6WFoPX694icaiYiIiIiIiIiIqNx4opGIiIiIiIiIiIjKjdFpIiIiIqrz9hRlIhgIKFFpADhQFI8O5hWLqGCRGhUMF8ffNttFalTQCIu4tMMyoMalDZeooIxSG1qmynCKCuqxQSVuJz5boGW3lNif6RarFvtwjUS7xaNlrFobdIxEa5FCp5iyW1TQNQLosD/oEW6XqKBbrFqJS7tEBeW6Szxarvu0SJuMSMuxgBYDThPzgkqMum7HhcvC548cfVIttj9SH6FIALmRdGX7wXC8HuaG4zHqQ6I2AkCBaCUhayMA2LKVRFgs65Foj/VQRqcT6mHUoQbq9dBrrFrWGz3PajrUVD0SLWuPz3me7dZKwiHqrEesHaPOejpWrru0krBd5jnWW9d6qI851ED9eB3i0no9lHFpGYEuGUseidbrpileVH1iHwnRabEPUxtzikvr87xyizdLbjFqp3243yZxm5ni6wI/0UhERERERERERETlxhONREREREREREREVG6MThMRERFRnfdTcQMEitKUqDSgRgIPF8aXiwoDyryojEvrkWh5pelwfLsRhjpPRgUdYoOAFhXUx5SoILwRMbyEK6TKOKDLVaeVCLPb1VNdY89ex5B0uSzzPF891WNU0PVK2D7nWLXjlVX1q7E6xAPlVaYB5ytL6+syBp2mXXVaxqWDvvgTMVjH48JlEa7jj8XeoiykBQLIj6otInKLQ2I5Hqs+XKRFp0VriUixlo8VcWlZG33FWt0UD7EZFs95l5qn1D9AjUTLH4eyXMbXiUs9VOum6TzP5crSjvUVcIw6u9Uy71eM9jaWOC95nXOthwktIsSyvJ1W5+BQAxOi06ash1rUWdZKM/lyyXry6HTCPKNuRKfLetXpSIq1kJ9oJCIiIiIiIiIionLjiUYiIiIiIiIiIiIqN55oJCIiIiIiIiIionJjj0YiIiIiqvP2F2bA7wsityikbM8XvcZkX8ZIkdZ3rFj0HStS/1Zvip6NyrLW2kj2ZTQjHnuS6b3GxKphp9D/yedzXDf0XmM+sS7GbG2e7RO9xuRtfHofxuTLJbfzNk/pQ+awfNT7KmdPMvdejur3xLEvozbPEP3ETNmfTO8TJub5td6LTn0ZZU/GkvXkfRlDPq2p6DHMV8cfiz1F9RHwp+FIRO3ReCQcXz8kaqPsyQgAxYXx0wh2oVpTZH2UPWvNYvUYlDGH2qiv6z0a1Vopflb02ui1Vup9FJUxU0xz6dEo1+W8hHro3PdWqanKPPWubIe+jJ770urrbrXXYR+23m/WrVb6HGqg3nvR71APtb60sgbqPWv9onb6Hfo1AkDAoS+jX5tnihdfvSej6dK/0WmeG7c+iin1XnTbn/6k0kQCqdVCfqKRiIiIiIiIiIiIyo0nGomIiIiIiIiIiKjcGJ0mIiIiojovtzAdPl9QiUoDelxavDXWotOGjAMWq1EjGQl0igPq60ZULmsxLJl0VdNbapRaWYYjGfPTY89wiwAaDtG+hDhz8mifrccQXSOAye/LLRKtxvzKENN23Ic2zynOnRCPFitaMl15rMTtjITodPK4tE+bJ+OBAS0S7RSXllFpwDkuzeh0XF2PTu8vqAe/GURhRD0dUFAcr4dFxfGxcJE6zyqIrxt6K4nC+M+iTyyb2kPqvW7KNhN6dNohLu01Kq3VDVgun8OSc2ULioTaI+c5t5yQt7P1WLVjJNr5vtR4tPZ1OcSj9TH1vrQ6p9TA5Msl66LlhF/7Psh6JuPRfi2mLOb5RF3za/Nk1DlNayXhFzVQxqP1uuk3ko/J7QBgGi6xaodItIxbp8ot2uw1Yu0+z/2zh2E/o9NERERERERERERUTXiikYiIiIiIiIiIiMqN0WkiIiIiqvMOFQbhM4MoLgoo25WrS4tlPQ7oKxIRQP3qqcXJ44EJUcGwiAAqV1lV58l4tH7Vaa9XmpaROhmdNnza5wyUCKB+NWl5BenkV6DW78tp+ehjTsseI9EuV0j1HBXUYspOserEq7GKqKC+D4crqxraFVJNU0YFvV1ZOuBzvnqq05WlAee4dHodjwuXhU//4a1jDhzJgA9BRCJq7jUcjq9bxfFl/crSsj7KeDQA+GV0WtRKX5F6DL5iUQ/DyZcBwJBXndauSK3UymgK0WmdHqWW9+V0NWmXq04rcemEFg5lr5WubSBc4tGer0jtdsVoZcy55inr+hWpZYsIhytLA2pcOhCQNU+PR4votDYm615A9CPR56WJF2DlStWGXl/jt9Mj0epVp537mPg8XnU6mkLUWY9HO0Wuy7rvcKA4ycyj4ycaiYiIiIiIiIiIqNx4opGIiIiIiIiIiIjKjScaiYiIiIiIiIiIqNwqtUfju+++iwcffBAbN27Erl27sHLlSgwfPtz1Nu+88w6mTZuGL774Am3atMHtt9+OsWPHVuZhEhEREVEdV1QQgGmkwSpSe42hOP53d7NILqt9jGTLNl+x85i6rPVxUvoyip5OUb0nmVjWxiDaP3ls96T2EDPUzxkYSq8xrXeTHBPLtqHNk7fzeesnZvtcxlLoJ+a571jCcdjO85SeZx57kunfFKUvY/JlQO1J5hN9x2TPMKAMPcnEWEjrvejUl7G+3kTvGFZUx/tVFhxOg2kFYUW0J304vm6IZVOrebJnrd6j0RRPI1+hvI32nJf9Gx36NQJqrTSiWv87r30ZZa9bUQ9de94m1LnkPWvthL63cp5zL0e33oue+8M69Zt1q3mGy5jpUg+d6rdbf1ytzhkOfRl9ei9ah76MwYDabzYox3zqmFMNTDO1uiluJ/syBrXmybL3YsCIamN20mWdD8n7N0bL8Pk/vRej03anfSb0cnTo+Viq2J9aLazUTzTm5+ejR48eWLhwoaf527ZtwwUXXIBBgwZh8+bNmDp1Kq655hr84x//qMzDJCIiIiIiIiIionKq1BONQ4cOxT333IOLL77Y0/zHH38cOTk5ePjhh9G1a1dMmTIFl156KR555JHKPEwioiq3aNEinHzyycjMzERmZib69euHNWvWxMYLCwsxefJkNG7cGPXr18eIESOwZ8+eajxiIqKKx1pIRMRaSER1S6VGp8tqw4YNGDx4sLJtyJAhmDp1quNtioqKUFQU/3x2Xl5epR0fEVFFad26NebOnYtOnTrBtm08++yzuOiii7Bp0yaceOKJ+O1vf4vVq1dj+fLlyMrKwpQpU3DJJZfg/fffr+5DJyKqMFVZC60CPwC/EgcEAENEAGU80NTSQjIqmDAmIoCmXNbmqdFpcQxqCguGiPkZlvOYjAp6jgDq8WifyMBp0T65bvuSx6MBLQLoGsVG8nn6PjxHoj3GtF2jgmKgMqKCIoooo4KmFomWcWmfGJNRaUCNSwf06LSIAMqooB6dTneITqfX8bhwWZhV/FhU9ftCKz8AWAEYUfWHwwjLOudc82RrCT1x71Oi07bzPBGX9oVldFp9zpuR+LoRUccMS6yLWLVSJ93o8Wjljp3roeOytm471VBAbTPhVlNTqGUJ6Vp5iG5jbjFtJVbt3ErCqV0EoNVA2SLCr35fneLSQa3mhfyyzkUcx9KV6LQWvxbrclmPRwdE5DohOg3xXNZftOU8h7GjxZelqJhrJXwzyzcv8bgMFKUYna5RJxp3796N5s2bK9uaN2+OvLw8FBQUID09PeE2c+bMwaxZs6rqEImIKsSwYcOU9XvvvReLFi3CBx98gNatW+OZZ57BsmXLcPbZZwMAlixZgq5du+KDDz7AaaedlnSf/MMLEdU2rIVERJVTCwHWQyKqHrX+qtMzZsxAbm5u7N8PP/xQ3YdERFQm0WgUf/rTn5Cfn49+/fph48aNCIfDyie8u3TpgrZt22LDhg2O+5kzZw6ysrJi/9q0aVMVh09EVCFYC4mIKq4WAqyHRFQ9atSJxuzs7IReE3v27EFmZmbSTzMCQDAYjPWyKP1HRFQbbNmyBfXr10cwGMTEiROxcuVKdOvWDbt370ZaWhoaNmyozG/evDl2797tuD/+4YWIaiPWQiKiiq+FAOshEVWPGhWd7tevH/7+978r29atW4d+/fpV0xEREVWezp07Y/PmzcjNzcWKFSswZswYrF+/PuX9BYNBBIPBCjxCIqLKV1W10CjwwYBP6cMIaH0TZY/GYmUaZMu2hH5lxXKe7TxPjkXFckTr6xeVPRq1/ld28h6NrpS+iT7nMb1fmZl8zDbd+ol5nacfo5hneOy96LmXo3ZfHvuayd5jrj3JxFhCj0bRh0z2ZfRp82SPRr/oQ5bmV/uJKb3LtJ5ksveY7FcWNNV+YrJfWYZ48mboT/pjmKn/8FaBiq6FgHM99B/2wYz4EnvAynro1qNRrCf2aEzel1FuB9QejaZcDms/Q6JnoxHVDjjqUA8996x1+dyVPuaT/RbNpNtLxpL3s7X1eW610udxnkPP2oR+s7K+JtRKO+lYwjz50qH0rHXp0aj1XjRF/fL7nfvNpvllnYsvh7R+gRn+eM3SezTWE2Oy92K6T61zciwknthB7Ukv+zL6YGtj6n2XMg2Pr9EurISGm3FR8Y3W+zxGHfoyep1XqjCQ/Gs7mko90Xj48GFs3bo1tr5t2zZs3rwZjRo1Qtu2bTFjxgzs3LkTzz33HABg4sSJePTRR/G73/0Ov/nNb/DWW2/hz3/+M1avXl2Zh0lEVC3S0tLQsWNHAEDv3r3x8ccfY/78+bjssstQXFyMgwcPKn+93rNnD7Kzs6vrcImIKgVrIRERayER1R2VGp3+5JNPcMopp+CUU04BAEybNg2nnHIKZs6cCQDYtWsXvv/++9j8nJwcrF69GuvWrUOPHj3w8MMP4+mnn8aQIUMq8zDp/7d3/0FylPedxz/dM/tDQqwECmiFkYjuLIMAC/PDxgu+MsbYCjgUxCqujsIFjklSdgQBdGcnuooTjG0EdhyME1kYm0POOTKcHMQFVwDLsiWXyxIGARWwzwqOCcJGK+ou1uoH2tnd6b4/1ux8n6fneRjN7uxIu+9Xaivd8zzd09sz81156M+3ARwRsixTpVLRueeeq46ODm3evHlsbOfOndq1axdXeAOY8qiFAEAtBHD0aukVjRdddJHyyGXL69atq7vNM88808KjAoD2W7VqlS699FItXLhQ+/fv1/r167VlyxY9/vjjmj17tq6//nqtXLlSxx9/vHp6enTjjTeqr68vemdBADjaTGYtLA8mSpUoKUQA68cDY1FBfywUly55EUAbkbbLNirtrxfGbJTaSxG6EwNxKC9650SivWifjQ7mNkYY2YdiEb1SZCwU2fPjzIH4dWx//qUVmY0tNxixdqKC/qm1UcHUj06b17xko9Pui1c2sWobG+xI3Xk2Lt3pRaJtdNDGo/2ooI1IzzLZ1pl+BnYaS0vNxQWbNdn/LiwdTFQaSYrRafOWctpKxKLTQ16dG6o/Vvaj0xXTVmDILA+7B5WM2Oi0Xw/N3Fhc2tYsuw+/Tjqxar+VRPrGy956LPaclSJR5ybqXLOtJILtI7xa5kasG20X4beIMO8HM6+j7NUy0zLCxqVtVHp03dY59016TNnUNtsiohCdrm3XnYzUfVySOs2YjVFLUsn8MU79D5Uzr7EodSzObKPPNjpdjUSsM2de+FrDap3rEDubrIVHVI9GAJguXn31VV177bXavXu3Zs+eraVLl+rxxx/X+973PknSXXfdpTRNtXz5clUqFS1btkxf/vKX23zUADCxqIUAQC0EMLXwRSMAtMF9990XHe/u7taaNWu0Zs2aSToiAJh81EIAoBYCmFr4ohEAAABTXnooUSlLnDigFI4H+jeRTGN3kzZznRh19LnC8eg0dtdps540etdpE8NLIndS9WOEeVI/2le8Q2pSd7kQe47eufqNlwv7jMSjYxFuTXBU0G6XeGPunabNsheJdmKEZqxwZ+nIHVhDcWk/UmijgzYufWw6KIwq+x/eKabjoFQaUTQ67dyBesR7X8daRARi1aUh7w7EJjpdsneWHvGi0/ZO05l3wI3WQGeHkbvsOndu9u8mXf9O007Nk5SX0/rL/rxSLFZt5tl6VXKmjb9u+nMbbSVh66E3L9QuQnJrYNnUsk7vrtP2LtS2JYR/Z2lb22xUWvLaQqT1a54kdZt+KjPT2linF4+2d6ROvb4ldq4/ZpUCsepYnDmL3FbFRqz9fYS28+e9UUy77LXoaFRLbwYDAAAAAAAAYHrgi0YAAAAAAAAA48YXjQAAAAAAAADGjR6NAAAAmPLKh6RSVuybGOpDFp/njrnbRXo0On0Z628/+lxm3evR6PRszAPLviTcd8wZ83uXlQLb+fOcfoiN9WEs9G90tgtsExlrdF50H17/M4WO1+vl6PQk88ZKZizUn0xye5KFliWp2+nDGO7R2GXefLY/mb9ul49J3d5l01nSZF+yo0X5tVylkTzaozGN1Kg01r/R1sAhu5wF5yWmR2Na6NEY7merUD302ZpVivRotH1ky15BSG2/xVLdZSlcA2M9GrOyXysD+yjUslAfXX+eXfZ7zNYfK+6jfl/GxOvDaGtgye/RaPrPdti+tF6ds71pbS/aYs2rX8v89WNLg+Zxr0ej6b1oa2CH16jZ9nIsyT2Hdq4/ZqWBHo2+LNKz0enLaF4If5tQ78VGezm+vo8SPRoBAAAAAAAAtAtfNAIAAAAAAAAYN6LTAAAAmPLKh6RStRi9sxFAN87sbp9WI7Hqaig6HY5Eu/O8qGBmo4Je1KrRuLSzw3A8Oi+Z6w5K4Vh1HohH++vOPC95GI0zO9Hk2HPV3yYvRLHtsnue8kAMujDPxJ5llm1UWpJSG532o4I2Lh2IDUpShx0zUbVO783WGRnrMhHAmTZSWHKjgjY6aKOCfqRwWkumdnS646BUGpYSr4bYZKcbo/bm2bFhv6aaz0MsEm3GEjvmzVNm1v2aF6qBfnuHkFg99NtM2Ci1jUGX3Xm5WbfL8Xi0vDF7TPUfHz3GwDaH00oiMM9vEaHAWOK3kjC1zG8lYVtG2HrY6cVzbW2zyzZSLblRalvzRtdr9awrEI+W3Oj0zKQ21unVABud7vDGbFw6TfK6jx8W81bxI9BZntQf897ydszGqqvee776Btcepmljce/Cdk1tBQAAAAAAAAAGXzQCAAAAAAAAGDei0wAAAJjyXr/Lqn8DRTemXHu8MM/Go/0Y4Uij8+xYVvfx0X2YMf+u0/5dV0NCcelIPLoYI7R3TLWRQgXnOVE+f3+NRgUDMerCeiQOGJoX3c6fF7rTtB8VNFE5PypoY2clM9YZu8tqYFlyo4Kxu06H7iw9ul4/Os1dp40pftfpUiVTOcuKd502b1/3bs/evEbrnI1He20gbFzarXnhu043emfp4k17A1FqP2Jt49LeXadDkeji3aTrjxXaQJhvYor7MMvNRKL9X7eZWlmYV78GJt48WwP92K2tgW4rCfcN5owlNmLt3QnaxJ5tPFpyo87dyVDdbQrz0uG6j0tSp+ofkySlJiJdmojotOFHp0N3k84i8zJz7KHtQ2PVBu+U7eOKRgAAAAAAAADjxheNAAAAAAAAAMaNLxoBAAAAAAAAjBs9GgEAADDllQ/lKo/kxZ6Hpv1QGutJFujlWByrv1wca7APo9+T0fYoywKPe2xvsCQL9yTz+4Q5/cvSSK8xsxrrw6jAPH9uUz3JCvPy+vMiY7nXXzHUkywteX0YS6YPY8ntZ2V7jZVNHzK/J5ld7zR9Gbu8Ho1dpkdZl9evLNSX0e+9eIydl5gejYnby3Faa7Iv2dFirB76dcP2aMzscqyWNVjnRrzei4G6Ge3D6PN7LB7uNn7P2kg9dPsy1l8enWd645llf3+ZWc/cdpDBOtp8L0dby8Jjbi33Xgfbztf2aEy8emjGSt5YyelZW1tOvXll8+Yrm36pHd4fX9sb0R+z652J3ceIN6+23imzjTJvXrhHo+3FaE+v/ye1Gan3efDP6ev8/oqpeTEz855PIz0aS3lxbEj0aAQAAAAAAADQJnzRCAAAAAAAAGDciE4DAABgyut4LVO5IytEABWIBxbigHn9eYV1Gxv0n8vuIxKddvfhxZZMjKoQe2xE6kcFbR7OjzOb+J6N1MXmNRp79hNaznb1nze2j8N5LjdyHY4UykSkk1KDUcHUi06bKHVHapfdmF+niUF3mrFOLx49o1Q/Hi1J3emwGavUXZbcKLVdPjYlOv26JJ3a0el0JFOqTIUUpq2BtjbmjdU8fyy0XNiu0VrmR6XtZzaLzQvUOb+W2Xi0VyuduLSZl3W482wk2kadc++bF2fMqz02Sh2KR/vrzdbeYD0s1GjzWkai03a9UCvNuhOP9mLPNi5to8J+xLpk9mHjy/5YGmmF4G8XYmPa/jahuHSzV/XZoy3ErwPdU3xOxDryK9rItX9+szpR6kZxRSMAAAAAAACAceOLRgAAAAAAAADjxheNAAAAAAAAAMatpT0af/CDH+jzn/+8duzYod27d2vjxo268sorg/O3bNmi97znPYXHd+/erd7e3lYeKgAAAKaw8mtVlcvVYn/FPNDwKNaTzOt35O4j1qOx/nMVei3a9UKvyPH1NfN7d7ljkf6NkX5icnolxno+1l+WvL5haWA5to9o3zHveENjfsO6QE+ytOR2xiqZdb9HY0epWne5s+T1aDTrXaYv44zSsDPPjs0sub0XbS9G26/xmDQ8z/ZlnOn1SZvOqpGeblNBOpwrrVM/3PpiHq+Ga1S0Z61fA0NsrfCOy/aELXSMSw6/IZ7tjRjt0VhKw2Mddl5sH7WxzJvn9mH09lGy85K6j0uRGhipm8Xei4Flvx4m9ZdtbZT8FpjemA5f2mAPxaq392rhD1Vj24VkZl5xm8C/D5psceh/3NzjMPMiT1A1f1Tjx272PY6ejL6WXtF48OBBnXXWWVqzZs1hbbdz507t3r177OfEE09s0RECQPvdcccdSpJEN99889hjg4ODWrFihebOnatZs2Zp+fLl2rNnTxuPEgBaj3oIANRCAEe3ln7ReOmll+ozn/mMfu/3fu+wtjvxxBPV29s79pP6/3XVqFQq2rdvn/MDAEeLJ598Ul/5yle0dOlS5/FbbrlFjzzyiDZs2KCtW7fqlVde0Qc/+ME2HSUAtB71EACohQCOfi2NTjfrbW97myqVis4880zdeuutuvDCC4NzV69erU996lOFxz9+7UdULne38jCPCqXFB9t9CEeUfaeOvPGkaeKJt9/b7kM4Yuzfn+ktbXjeAwcO6JprrtFXv/pVfeYznxl7fGBgQPfdd5/Wr1+viy++WJJ0//33a8mSJdq+fbve+c53HtbzUA9HUQ9dPadTWwAAIABJREFU1MMa6uGodtVCaXLqYfnAkMrltJAGC8aP/ccjScqG48yhmPaE7M+PuTUWgXLi0oW4nYns2Yi1H/Nz4tL2cW9/kQh3KBIdi2nLmefHnsNjdt0ZK3kxPxuXDixLbly67MWqO8xYR2qi06lbgzsDcekub95ME3XuTvyxWiTaxqWPSYaceXb9GLOPY9KJi80d7bI2nYvJ+rfh64ptGyJjDXKizvbz6+/ORodNaj/xLzIyRTsv9EgIHKNf/5L6tUyRWmYj0JIbpc7sctmLRJv1rCMJznPj0QqOKVIP81Jg+XDqoa3Z9g+kfwptK4nIx8PGpZv9GGXmF7WxXz/aa+PRmXdyquaXHs7LZtltEVEyb8zBSEy7qqpZdo/D7sNGvUsNxr5jqpE4cywSbc9VNGIdjVKnquTNtZE4om4GM3/+fN1zzz36h3/4B/3DP/yDFixYoIsuukhPP/10cJtVq1ZpYGBg7Ofll1+exCMGgOatWLFCH/jAB3TJJZc4j+/YsUPDw8PO46eddpoWLlyobdu2BffHFd4AjlYTWQ+phQCOVvzbEMBUcERd0Xjqqafq1FNPHVu/4IIL9K//+q+666679D//5/+su01XV5e6urom6xABYEI88MADevrpp/Xkk08Wxvr7+9XZ2ak5c+Y4j8+bN0/9/f3BfYau8AaAI9lE10NqIYCjEf82BDBVHFFfNNbzjne8Qz/84Q/bfRgAMGFefvll3XTTTdq0aZO6uycu0rxq1SqtXLlybH3fvn1asGDBhO0fACZaK+phqBaOzOqUyp3xOGAsNhi7I3UeHmto/4U7YUf2F9q/l24Kxq/9O7p21v7nQNbp5veyrtp6tdNE1DrcUFTWaSJaTlTQO8QO87zeWG5ihaE4oOTFqmPzIhHu0N1TY2IviRNRy9wdVqq1A3MjhY1F6rJIHK7incRBs/5aVrsQ42DJvSjj2PTQ2PL+dNA8PiiMOjAyuXednvR/GyaSEi9GLHlvdLPsx2htSfH2Ye86ncfm2TtX23h0o+0iYqLRabM7787ScqLT/l2ibQ0Mz7MfyyxSo2J3k86cu1WH9xFsOXE49dBGqdPA45L7oif1X7vD4Uaivb8p5oCHzQkYTt1fbNj8ooP2D4yktME7x1dNvr9qXn8/UtxhotOdXvzaPlcpdAdqfyx2TA1Gnf3zFppnY+SNbvO6g97v2qgjKjpdz7PPPqv58+e3+zAAYMLs2LFDr776qs455xyVy2WVy2Vt3bpVX/rSl1QulzVv3jwNDQ1p7969znZ79uxRb29vcL9dXV3q6elxfgDgSNaKekgtBHC04d+GAKaSll7ReODAAf385z8fW3/xxRf17LPP6vjjj9fChQu1atUq/epXv9Lf/d3fSZK++MUvatGiRTrjjDM0ODior33ta/re976n73znO608TACYVO9973v13HPPOY/9/u//vk477TT96Z/+qRYsWKCOjg5t3rxZy5cvlyTt3LlTu3btUl9fXzsOGQBagnoIANRCAFNLS79ofOqpp/Se97xnbP31y7avu+46rVu3Trt379auXbvGxoeGhvRf/+t/1a9+9SvNnDlTS5cu1Xe/+11nHwBwtDv22GN15plnOo8dc8wxmjt37tjj119/vVauXKnjjz9ePT09uvHGG9XX19f0XQUB4EhEPQQAaiGAqaWlXzRedNFF0R4L69atc9Y/8YlP6BOf+EQrDwkAjgp33XWX0jTV8uXLValUtGzZMn35y19u92EBwKSbqHo4PKukvKOkQtumUN9E/5+w9dtT/WYsrzuWVP3+io09l9u/sbHnKvZyDO3PnWf7MhZ6NNp+i6Yvo+3JKHl9GTvqPy65PRszv/9ZoK9Z5vcas/3KTA8xv+2Us17oSWaWnRdMYaZnWOb1YbTrI1X3gE3LLw012KPR9icb8U5ApVQ7UZWy16PRnPzX0lpfxv2Z2/NvZjpzbNn2azwmHQoe03Tz2nBV0i/bfRiOify3YVZORz/Tkf+tntjGfoU+srbpn7ddZmtUuJbZ7wmSkjPgTQweYuOcHoVmJfH7K9reiP5YrXDYOpR1ROphuf7y6Hr9ZX//tudjvGftBNRDq8H+tTFZ5KV06px3wEOm7pXNclp1+zCWCn8gGzkm9xezfR+dno9JpzOvIxmpPa/3Dwnbv9E5vgb7RMZUIz0VM/PiRfs6mn1kke6J9Z7rtWpzPRqP+JvBAMB0sGXLFme9u7tba9as0Zo1a9pzQADQJtRDAKAWAjh6HfE3gwEAAAAAAABw5OOKRgAAAEx5I8ekyjvSOjE/sxyJIjsJ2wbj14mfOGo0Om33HzkOm0vzjym0fz/OnXWa6FWnew2CjRFWO208OhwBtLHB3PtfGm6M0B1rOCpYCiz7qTEnVu3+znmjcemQ3I9Om3Po7W/EOd/h/+ll43xVs78R7wTYSGHFO4ldaS3aN6M0bJa7gvMGTIy6Ox0WRg0Oj7zxpKNY1pH85nPsvWGDNSo2z9t5qL2DH792am8kHx06pojci0Q70WkzVogYm7h07kenS6FItLsPZ6wjUvNsnSv7MW37vGa5cLyNzYvWQ7ueRmpjg7UyN++V3Hvf2No2bNpMlL0/YLa1ROq3IAmoetfQ2Rh0xZzQrsT9bNu612H+aHd48+yYH4lOTYS7ZI69mWi3z/+9nDHb0uMwI9GN7D/LEx2qNlcLuaIRAAAAAAAAwLjxRSMAAAAAAACAcSM6DQAAgClveGYyerfkQsTYLDuRZTfyFYv5hbYrxpntNo1FpwtRQScGHZkXiIQn/l2nbcyv049Emzhv4E6qo+tmuUPBeW5U0DtcO+bcSdWbF7izauxurIXIX2gsEg10UqSF6LR5zb3I5khS/7qOzDsBTqTQ3AV1yLvrdGda226wVPXGahG3g6VO87gfna6a5RGzTHT6dZXK1D4XI92J1FF8w7vtI5L6j0vxOHOgBvq1x6l7DUaxm2bj0pHPvBuddseyQKw6Gp0O1MbYvMJ2pfqPj66bethgxLpwqVnofHgvRFJ4YX6zb68e2pfL1rXR9dpgNa29OSrVxtpK+HeMHslsrXT30enUttrJ7/B6mtga2GFqY+r9YY5Gp816aUJukx4Wurt0FotHR+9cHc/ED440Vwu5ohEAAAAAAADAuPFFIwAAAAAAAIBx44tGAAAAAAAAAONGj0YAAABMecMzEmVdSaH/l9NqqcE+jEnmN/YKzVNknu1j5j1XYF5hn+YYY/Ps75JUvf6CjfZebLSfmJlX6MPo9G90x/JAH7K8lAfnOW2nUm+eM6awJNYczu6w7uLoujn3Ve/JbP8y29ctLfR5tD0aay/eiNfjbMj0byyNuEdSMtt1ml5j9nFJKjvzRuo+Pt0NV4bafQgtVe1Opc60+GY2nBro181IjQr1W4z1aHR7QzZ+HE0xh5unXn/BNDxmeyU69aoUrpt5oDZK8XoY2kehF22gHhbqpqmPeRp7Mc3finjrvqA81lOxautZ7RfLvH6zdrsR84v59dCpZdWR4FjZ/EEsp+5zdQTGSt7fg7Lp0ZhG/laU/EbQE8z/G/M6/1yHxHo5FuYq0VCTtZArGgEAAAAAAACMG180AgAAAAAAABg3otMAAACY8kZmSnlXMR0bjk5HIsteMioUly5EokOx6shz+fHr8HPF5tmsoB9FjkWi648VY8/15xWi041GAAPLUjgeWEiDxaKCdr2Jyy5yLypq03x+ND3PM7Ncm5gk7hNXTcwvTU2MuurFo51oox+dzs1yFp5n1suRedPZyGuVdh9CSw3PlLLOOt0CnBRtUvfx4jx3yK1t9nGvzjmtH+pv80bH0cxb1im3fuo7rb8suRHpWI1y61ykvjbTSqIcaSVh66F3TLLrhVoZWG7w5Hp/UpTZtiBe/nokGPt116u2BpoCW0rcE2XrVyntcMZsvLnkxKOz4LzUvMHSJDIvcm5iYxOt8bh08/OGB4lOAwAAAAAAAGgTvmgEAAAAAAAAMG580QgAAAAAAABg3OjRCAAAgCmvOjNX3p3X6WVoehLFeh5WzXKjfR79nmQN9ldUw30eTS+sqiLzws9l+zDGe43V36Y4Zh73+46l9ef5+7B9yKK9HAP9GiV5vcbcIWfd6XmooNz0rkq8RnFOz0a/55vpL5bb18F7E1UT29csPC8+Vv93ic4zj09mb7EjXfXQ1L4eZ2Rmorwrifc8DC3Lqyn+2yZQvwo1NQ/U3gZ7PtZ97sPkt6SL9Wi0NSVzeiM2Vw+dfcT62Zp6mBXqYaBWevXQ6VMb61kbO6Gh+uidRNuzsVpNvTFzvLamZu68xOmvGO5Lmzh9ExUZC+8jDdRD35HSlzGk0T6Mh7uPZvvVTu0KCgAAAAAAAGBS8EUjAAAAAAAAgHEjOg0AAIApb2RmrrQ7V1L148wmNmUTZF4UueFIdLX+45Kc2F9q91+Ic0eeK7D/aMTa2SYS8/NjeSYi7UYF/Xn19+FHD0Px6MJYdB+BqKA/r8GoYBKIURfkdtHPW9ZdLOzfe4UUWw2JxbvdiY1F+cYftpuastfafQStNXyMlHXVeZvYGthgZDna+sGJx0bmxdpWZOEPWDDq7bHzoglT+3kt1JT6bSYKLSdCtazB1hRSOC5dqJvOPvK6y5LcGujXQ6fNQv3HC8xJ9F/WNKs7bXSubfcRae/gHF7kOGLbjXebhmttk/tvVj4BEenD0WwbCa5oBAAAAAAAADBufNEIAAAAAAAAYNxaGp1evXq1HnroIf3sZz/TjBkzdMEFF+jOO+/UqaeeGt1uw4YN+uQnP6l/+7d/0+LFi3XnnXfqsssua+WhAgAAYArLZmTSjKwYU7ZR6kbv4uyPhSLX3nOl5rmyZu9wHYpOFyLhoW38yHL9ePTo2BsvS/7dn8PzslDsObqPcMTa3lk1j8QBG73rdMOyxqNrDd8JtOF5DT5xw/MIT9eTHfJ7J0wt1Zm58u68GEW2Kw3fWXr8rSTi0elI/DoW4Q7Mc3ce2aQQnQ4sN1sPy5F6GIhLxyLWTq3066Ez5g4lgVh1wxFgr4ZkkbFG9xGcFjukWOcL6lzTskPNRcFbekXj1q1btWLFCm3fvl2bNm3S8PCw3v/+9+vgwYPBbX70ox/p6quv1vXXX69nnnlGV155pa688ko9//zzrTxUAAAAAAAAAOPQ0i8aH3vsMX34wx/WGWecobPOOkvr1q3Trl27tGPHjuA2d999t37nd35HH//4x7VkyRJ9+tOf1jnnnKO//du/beWhAsCkuvXWW5UkifNz2mmnjY0PDg5qxYoVmjt3rmbNmqXly5drz549bTxiAJh41EIAoBYCmFomtUfjwMCAJOn4448Pztm2bZsuueQS57Fly5Zp27ZtdedXKhXt27fP+QGAo8EZZ5yh3bt3j/388Ic/HBu75ZZb9Mgjj2jDhg3aunWrXnnlFX3wgx9s49ECQGtQCwGAWghg6mhpj0YryzLdfPPNuvDCC3XmmWcG5/X392vevHnOY/PmzVN/f3/d+atXr9anPvWpCT1WAJgM5XJZvb29hccHBgZ03333af369br44oslSffff7+WLFmi7du3653vfOdkHyoAtMxk1cJk5oiSGSPKvV6Gue0vZsZyvw9frB+i0+cxPC+rmv5Xob6O3naN9oos9jULzPN//0Z7L8Z6kpVC88L9FaP9zyK9HG0PMacvo3/5RKRfme1JZpcLfbzsZmas0CfMblcYG+c8Xx5YlpSE9t/gczXTrnLKGvTfeK03mf8uHJlVVTqj6r5nJLfBntP/MFYP3aFQv8VoL9pIf1xF9hF83th7vkF+jbL9HJ2PWrP9Zs03MYV6aOZm9hubQk0N1Dlvnq2Pid+/MdSXMVaGYucz0sM2D9Yo7++SHctidTMyFjqOBt8bhc9GTLtq52S0nmyyX+2kXdG4YsUKPf/883rggQcmdL+rVq3SwMDA2M/LL788ofsHgFZ54YUXdNJJJ+k//If/oGuuuUa7du2SJO3YsUPDw8PO1d2nnXaaFi5cGLy6W+IKbwBHJ2ohAEx8LZSohwDaY1K+aLzhhhv07W9/W9///vd18sknR+f29vYW+k3s2bOn7n/dkaSuri719PQ4PwBwpDv//PO1bt06PfbYY1q7dq1efPFF/af/9J+0f/9+9ff3q7OzU3PmzHG2iV3dLY1e4T179uyxnwULFrT61wCAcaEWAkBraqFEPQTQHi2NTud5rhtvvFEbN27Uli1btGjRojfcpq+vT5s3b9bNN9889timTZvU19fXykMFgEl16aWXji0vXbpU559/vk455RT9r//1vzRjxoym9rlq1SqtXLlybH3fvn38gxLAEW0ya2HXjCGVZqaqVt3/zp6Z9cxErQoRazOvGL+2O7QRQC/XFIhEF6PYJsrWZNwwFKtOvXmheHRhLBZ7tlHnyP4Ui1WHtivMC8SlI3HA6Jg99ZH4Wx6J+Snw+hf2GYn5OdHUQHw1Ok9eXDSLRKJDEVOi02OSwcnIJNa0ohZK4XqYzBptJeG/5k7LCFsPI9Fp/z3faCsJtwaa5QbbQIyO1f9MFSLWzbzP/Y9yWn8sWssarZuFdhQ2Bh1+rmCd8+YloZqnxuPS7gGaxQZbTkj++8sMFN5fodf1MN6Hoci9Py/03vBrryKOhNo5ASWr3q+RDDZ3bWJLr2hcsWKFvvGNb2j9+vU69thj1d/fr/7+fh06dGhszrXXXqtVq1aNrd9000167LHH9IUvfEE/+9nPdOutt+qpp57SDTfc0MpDBYC2mjNnjt7ylrfo5z//uXp7ezU0NKS9e/c6c2JXd0tc4Q3g6EctBICJqYUS9RBAe7T0i8a1a9dqYGBAF110kebPnz/28+CDD47N2bVrl3bv3j22fsEFF2j9+vW69957ddZZZ+lb3/qWHn744egNZADgaHfgwAH967/+q+bPn69zzz1XHR0d2rx589j4zp07tWvXLq7uBjClUQsBgFoI4OjW8uj0G9myZUvhsauuukpXXXVVC44IAI4M/+2//TddfvnlOuWUU/TKK6/oL//yL1UqlXT11Vdr9uzZuv7667Vy5Uodf/zx6unp0Y033qi+vj7uOA1gSpnMWtgzo6LSTGloxM2ojWS1/+5uY9XxiLU/Vj9y7UesQ3e1TkbC0cPoHamzPDKv/jZ5JDrtX4KQO3d4No/7MT97+E4EOhwpLDyXc/fU8D6ceFj0Lqv17ywtuVHBxO7Pyxi7d0gNR5YViJv6c6Ox51AMPhYH9E+Ns13gef19RCKr01m1MrnR6cn+d+GMY4ZUmpkU7h5sa5ttJZF59dB+NrJCmwmz3mA9TALLo+v2ACNtJpqITsfutF646bAtB4EYtT+Wx9pF2Drq18pAPYzWudAdqBWPRycNvtXteyV492iF4/eSnNc82mYkELmP/z30x+rvL409V17/8cKY7yiKTh/OzbSl5mthS79oBADU98tf/lJXX321/t//+3864YQT9K53vUvbt2/XCSecIEm66667lKapli9frkqlomXLlunLX/5ym48aACYWtRAAqIUApha+aASANnjggQei493d3VqzZo3WrFkzSUcEAJOPWggA1EIAU0tLezQCAAAAAAAAmB64ohEAAABT3uzuQyp3ZxrO3AaDw9Xa+pBZHvZ6ko2YsarXo3FkxPQ1i/VyHLH9ympjeTnc46zQu8r2JLPzCn0ezbLp5ZiOuLsL9RMbnWzGQv3JJKcfWHx/4T5hTi9Gp/+Z35MssD+/J1msR6PZh+1dljfcrMzv8RXpvRjsQxee5/ZX9OfVlv1eY25fzsb2EVqe7pJKu4+gtY6b+ZrKx1SVee/lqll3a547z/a2HfH63tr+tk49LPRyNDXQqY3hz1exHtb/TBXe8/b4A/0a31BSf7nRXrR+j8b4WKC2lRqsc34fRr/XrR0zy86sBpv55YW+tPV7dPrr0b9fI6o75r+uaWDe6P7NstOjUQ3Ni/Zo9Pvjhk7vRPRujLwM0Zdogno2pk3WQq5oBAAAAAAAADBufNEIAAAAAAAAYNyITgMAAGDKO2HGAXXMGNJQ5v7zd3CkY2x5yMSqK1V3XihiPTqW1p1X9ePXJmLYTKRwdKf1I4ZJ2Y+e1Y8U5qm3v0A8enSsfiS68dhzZJ7/XMnhR6fdqGAsUhh+rsTZLpZRCyx76360040HxubVX05jcUAvAmgjgfa5/Li8u/+87uPTXTo0EZnHI9e8Y/ar45iKE4GWpJE80Eoi0nJi2NvH0Eip7rxYxNqJW4+4+7N1LvfGnHrotBzw47x5/bHIZ7mgwei0My8Wjw7VPMmJSCdNtIgolNdgttdj/gjk/jZ54Lz5fzicNhCRv0uBZcmtWelw/e39sUIkOlADE68eprYdSSw6bd9fsb8BsVMdGmsw5ixFos6NRqwP87marYVc0QgAAAAAAABg3PiiEQAAAAAAAMC48UUjAAAAAAAAgHGjRyMAAACmvHld+9XV3aFD1Q7n8UPlzrHliulDNujNs+t+j0bbz9Hpa+b3cgz2LnP/27/t31gte/3KzFy3d5nX57FkezSaZfeQGu811sS8QrOqWJ+oQP/GpPBcjc1LIseRJvUbaiXeMeWhFb/vWBbpSeb0TayNFfuJmbFYP7FY70VnLG9oXmk4PG86Gxluri/Z0WJe1351dncoy936UjE9bEfMmF83bZ3za6Wth5WR+rVRarweVm1v27LbOM/pZxvoXzu6oRnLA/0apcZ7NNpNovUw0ocx1nsx0KPRr1FJqN+s37M21r8v0Do395sBhvoy+ucsq/+3x18P1TxJSofMmO3DOOzNi47Vr22F5wr0qU0KPRoDvRzlne4GS0diTnwee4EKG5qnim3W8LzwYJ40Xwu5ohEAAAAAAADAuPFFIwAAAAAAAIBxIzoNAACAKe9NnXvV3VnWa1mn87hdf61aWz5UdefZ6KAfIxw0UUEbI7QRQkkaLtcigDZGODTiRQqdGKE7Vi2ZWLWNDZb8iLWNVdt5hxFnduYF8nWFfTQxT156Kw3EAb19xOLRzlBhH+EoYlNMjC6JRKJt3C7xou5OBHCk/uP+uh+rdmLQTqQw9+bVHysNeXnAaSwdrr7xpKPYm7p/re7usoYzt0ZVclPLMlPLvHnxethRd/nQiBexNrHqQbM8XPLqoYlLR2PVNh7tRXZtxNp+XvPDiU6HNNgGohCPTuvHo0fXzbKtV6n7GU2SwLzDqGs2Op3ZGuXlbYOnJnYO/Xpo1m37iEYj0aUhb96QHcvDY5EWEU6biWEbj/ZeE/u75P6Y2S72HgqNxaLtsdcyTcLzAm+C4rzI/tV8LeSKRgAAAAAAAADjxheNAAAAAAAAAMaN6DQAAACmvN7yXs3sKOlgITrdVVsumWV/nolS+2M2OnhwpLaPQS86HYpVV0ruPCdWXYgR1o9VV73otF3PbFzaix42HHUObaNIhDlyh9TiPhvYX+yQDiMeHYoYVgvHF7htZ+wuq4U7ldaWo5FoOzYUnmfjgY2O+ZFCG5FOK7Xl0jDR6dclI1P7FtzzO/ZqRkfZiUdL0nBu7iad18Zeq3Y5814r2TYT7j4Omfp4YKS23F1y59l9dJgxvx46d65O3Xo4UrKxahOjTt06Z2ugc3fqQuy3iV4KTbaBCMWj/fUkso/Uv+N1YJ5VuJu0qXNJpMy5+7Ar/qDZs3d+3btO1x5Ph71WEmbMxqXTivtUtraVYmM2Ou1HrE3dc+5UXfXqoY1He7Fq53d2sujhs2hfokbj0T7nbtWR6LSz/1iuvs5liM3WQq5oBAAAAAAAADBufNEIAAAAAAAAYNz4ohEAAAAAAADAuNGjEQAAAFPeyR3/rmM6UqcnoyTty7rHlu3YQW/e/pKZF+lXNss0ijrgzbO9zGy/Rru9FO/fOGj6lQ2ntR5SwyW/d1lqlk3vslKkZ1Shz2FoXnAXhX5lDeyu/nMH9hd67lhPstjxOr3QwtPifcdMKy/bg2x03SyPJHUfl9y+jE5PskJ/Rbvsj9l+ZbYnmdtrrDxYe/K0YpYHp3ZfwsNSHXrjOUexuel+HVMqFerGoOmvaHs0Hkzd82H71L6WunXuQLX2Puo0zfZs/1pJKpsPTmhZklLzGS2l7tiQ6cto+xXafo2SVDWfS9u/sdij0Sw22K8xVjedFnoN9mH0x9JYn8fI/i33d4l1X2ywR2WkZ61THyM9a51+jX49HK6/bHstSvF6WK7YXrTmPeTVQ1sfkxHbr9E7KKdHY6R/Y+68iTRu3h8wpy9jGp7n9mgMbPMGzyVJeXW4zsQ3xhWNAAAAAAAAAMaNLxoBAAAAAAAAjFtLo9OrV6/WQw89pJ/97GeaMWOGLrjgAt1555069dRTg9usW7dOv//7v+881tXVpcHBwVYeKgAAAKawN5UqOraUar8XAZyTHRpb3pfXon3RiHXJi1VXZ5ixWqRwhhc5stFpG6vuLrnzXhsx8cVShzPWYdZtrLoy4v6zfigNRArT5q4ziMalnXnNRcUa3X8zz9XsMQUVooL1lyU3HmhSpE4c0F+3cemSl+C1kehyxYsRDmZmnlkedCOA6WDtydJDteXkUEUYlWZT+1ycUDqoWaVUg7lfN2r/m9vWwO7EfcN2J53BsQ7zRu+o1t57JfmRaBNTtVFh70PkxIiTjuBYaHl03bSPyGpjWdWthzZinPsfdBsXbrBFRDPx6MJ2jcav7aH6x2Tm+ZFwmwKOncNgrNqPmDttJry5TnQ6PM9tOWHOU6Fu2lrpHq8Tl7b1sOLVw6HaejJUe+/aGLUkacRs50Wnk6pZb3F0OrGR6DSQ0/fWEyc6HZ5X77nTanMtNVp6RePWrVu1YsUKbd++XZs2bdLw8LDe//736+DBg9Htenp6tHv37rGfl156qZWHCQAAAAAAAGCcWvpF42OPPaYPf/jDOuOMM3TWWWdp3bp12rVrl3bs2BHdLkkS9fb2jv3MmzevlYcJAJPuV7/6lT70oQ+LGh9uAAAgAElEQVRp7ty5mjFjht761rfqqaeeGhvP81x/8Rd/ofnz52vGjBm65JJL9MILL7TxiAGgNaiHAEAtBDB1TOpdpwcGBiRJxx9/fHTegQMHdMoppyjLMp1zzjm6/fbbdcYZZ9SdW6lUVKnULm3ft2+fJGnjNx9Uz7G0oDxz+zXtPoQjyjGKX007nXz4bZe3+xCOGCPZkKR1k/Z8v/71r3XhhRfqPe95jx599FGdcMIJeuGFF3TccceNzfnc5z6nL33pS/r617+uRYsW6ZOf/KSWLVumn/70p+ru7o7svYh6OIp66KIe1lAPR012LZQmtx72lmepp5xqduZmUQ+YLNaxWW35YOq27TnWRKz3ZzPcMTN3v4lYz/Ri2vvT2tgME5f2707dldYiWoe82JK9i+traS2+WPbuxlo2d5OtpLV/8g+lk/rP/6NKISoYurOqn4Yzd1n1I4A2Lm3jgLEIoHOXVS8ebePSpUPhu0mXDtWe2MajJSk5VHtfJq/V3rv5oUPCqDyb3LtOT/a/DY9NhzUrTdWVuzHSSl6rG52mDnUkbh0qmTd6aQJaE2R5apaT4FjsTtCN3iW6auLSfmo0M5/lQuq10eh0IOpcTLY2Gp2uH4+W3Mh5TBY59sSJ2Da0uzinVnqx30CbCf+u006ttMsjXjzaiVV7Y8Mmmm/i0XZZkpKKaR9h7zQ97L7nExudHvEO2L5Zqt6YM63Rdh+RSLSdZ+8a788LxaojEWtfniZStblaOGn/0siyTDfffLMuvPBCnXnmmcF5p556qv7H//gfWrp0qQYGBvRXf/VXuuCCC/STn/xEJ598cmH+6tWr9alPfaqVhw4AE+rOO+/UggULdP/99489tmjRorHlPM/1xS9+UX/+53+uK664QpL0d3/3d5o3b54efvhh/Zf/8l8m/ZgBoBWohwBALQQwtUzaJS4rVqzQ888/rwceeCA6r6+vT9dee63e9ra36d3vfrceeughnXDCCfrKV75Sd/6qVas0MDAw9vPyyy+34vABYML84z/+o8477zxdddVVOvHEE3X22Wfrq1/96tj4iy++qP7+fl1yySVjj82ePVvnn3++tm3bFtxvpVLRvn37nB8AOJK1oh5SCwEcbfi3IYCpZFK+aLzhhhv07W9/W9///vfrXpUY09HRobPPPls///nP6453dXWpp6fH+QGAI9kvfvELrV27VosXL9bjjz+uj33sY/qTP/kTff3rX5ck9ff3S1KhP+28efPGxupZvXq1Zs+ePfazYMGC1v0SADABWlEPqYUAjjb82xDAVNLS6HSe57rxxhu1ceNGbdmyxbn8u1HValXPPfecLrvsshYcIQBMvizLdN555+n222+XJJ199tl6/vnndc899+i6665rer+rVq3SypUrx9b37dvHPygBHNFaUQ+phQCONpP9b8NZaaJj00SDudvns2R6yKW2n5zXe7SamOuVvEuXqqYvX1W296I7cTit9ZebURoy87wejaYzYabI2AT0b7S98Wy/xt/spcF92OX6vRb9eanXY9c+cxrZR6PsmY+dpwlh9l843ECvW7+3bah/Y6wHbrR/40htQ6cPo7du+zX6PRo1UlvPM/8DYfZZ9casPDJmpzmfr0hPRdsr0psX7POYeh/YSP/GRFKSeQ2FG9TSKxpXrFihb3zjG1q/fr2OPfZY9ff3q7+/X4dMo+Frr71Wq1atGlu/7bbb9J3vfEe/+MUv9PTTT+tDH/qQXnrpJf3BH/xBKw8VACbN/PnzdfrppzuPLVmyRLt27ZIk9fb2SpL27NnjzNmzZ8/YWD1c4Q3gaNOKekgtBHC04d+GAKaSln7RuHbtWg0MDOiiiy7S/Pnzx34efPDBsTm7du3S7t27x9Z//etf6w//8A+1ZMkSXXbZZdq3b59+9KMfFQovABytLrzwQu3cudN57F/+5V90yimnSBpt/t3b26vNmzePje/bt09PPPGE+vr6JvVYAaCVqIcAQC0EMLW0PDr9RrZs2eKs33XXXbrrrrtadEQA0H633HKLLrjgAt1+++36z//5P+vHP/6x7r33Xt17772SRi91v/nmm/WZz3xGixcv1qJFi/TJT35SJ510kq688so2Hz0ATBzqIQBMfi3sTFJ1JalSLxOdOtnWWiwz8yKV1bRSG8vca5eGTCR6OK993VBJOpx5XSb3OpzXtqmk7lcUZRNTLSdu7LWU1LYrmfhxyYsHl0xc1PmOouT+/v7v0gw3Ll3/8cI8bx+huHSz0WlH7kdsA/ufkOdyV5NQXNqfl+V153kvvzfPO782Lh2JTjsRaWfZjQznJjpdiEeb6HRe9fbv7KTBcxqJOjuRaBuXTtx5eRqIRHvz7D6SpE6sPhspPtaAln7RCAAoevvb366NGzdq1apVuu2227Ro0SJ98Ytf1DXXXDM25xOf+IQOHjyoP/qjP9LevXv1rne9S4899pi6u7vbeOQAMLGohwBALQQwtfBFIwC0we/+7u/qd3/3d4PjSZLotttu02233TaJRwUAk496CADUQgBTR0t7NAIAAAAAAACYHriiEQAAAACAaaI7Kas7SSW5/dcy07MxM43zql5zvCFzvVJH4u6j08y1Y12p2/NuMOsw8+w2bv+7jrQ2Vjb9H0fXa3OruVn2ei3a/o22lV9WdfdnexR6rfGaaq8X669ou+Glsf6N9pjqtNBrhH0t/efK8gZ3aqdFzoWze39e4Lm8l9zty2iX/V6OVbvsnUPbv9H0aFTmPllSrT+W+30YTY/GwlioR2PWZJ/LSN/E3L4xJ6CXox3L6/RozJvs0cgVjQAAAAAAAADGjS8aAQAAAAAAAIwb0WkAAAAAAKaJskoqK1WHl221Edshs1zyMqudJqbc4cWqbVzaxqhL8iPRtXlp1lk7ttTdX9nEWctextauD0fiwaH4sT8vHnRtLGIcij37W6exWHWDcWl/u9fljcahI/tohYafKhC/TrwocmLy7IUxu25z737s2Uapndiz+z504tJVf6xadyxvMjqdmNevsAcTfU5K5rpBL9vvxKBzs43/IuSRWLUk5VnxsQZwRSMAAAAAAACAceOLRgAAAAAAAADjRnQaAAAAAIBpopSkoz9eoDc1651mueLFMm3st+RFK0tOhDmru40/r2Tm+RHr2D7sWOwOz6GYciGy7MRIm4y9tjCKPJkx54Y1eUjRu1M389z+PkLRaf/24WZeLB7t3pHaGwvFpf3YcejW5f6dpe3uvThzkmb153n7Tkrmjur22P27U9tDTOvEpIlOAwAAAAAAAGgXvmgEAAAAAAAAMG580QgAAAAAAABg3OjRCAAAAADANJN61x35PRtr8+TNy+suj86t39PN773obtOe3oOt7nloz6bfXzLeU7L+vIaf198mr/+6ttpEPGvi9Ff0Bhs9Nfat5/dJDPVN9DU8zzzZhOzb/dzkWe3TmDg9FUvOPNuXUbYvY+btz4wl9T6ijf4OHq5oBAAAAAAAADBufNEIAAAAAAAAYNyITgMAAAAAgLbzI8YTsV2jEV67j2qb4satdqTEqo9qmXsO86w90f8jGVc0AgAAAAAAABg3vmgEAAAAAAAAMG5EpwEAAAAAmOaqDd7Gt2rCyFUvmJyN81qmrMkob2y70G+Ve9s0+9yNPK+/71KDEXF7jBNxl2z/d26llgeKJzP1nZgnS90nTsw5zasT/FyJ+3lK0vBYQ9LD2CZNmo7Wc0UjAAAAAAAAgHHji0YAAAAAAAAA48YXjQAAAAAAAADGraVfNK5du1ZLly5VT0+Penp61NfXp0cffTS6zYYNG3Taaaepu7tbb33rW/VP//RPrTxEAAAAAACmjWqeqZpn8v/PnZOP/WSS89OoLE/Hfqpyf5x5SsZ+YvvI8sT7qY0dKfI8GftpdJ7/k+Ua+xnvMUxmT8ZWyJNk7EeJ3B/LH0sT86PaT5KEf+w2Ser+2KdKEuenMHfsJ/Jc0Z/aPpI0cX6c/ZvjLRxTmtZ+7L790+ZsEzgHTWjpJ/Lkk0/WHXfcoR07duipp57SxRdfrCuuuEI/+clP6s7/0Y9+pKuvvlrXX3+9nnnmGV155ZW68sor9fzzz7fyMAFg0v32b/924Q9CkiRasWKFJGlwcFArVqzQ3LlzNWvWLC1fvlx79uxp81EDwMSiFgLAKOohgKmipV80Xn755brsssu0ePFiveUtb9FnP/tZzZo1S9u3b687/+6779bv/M7v6OMf/7iWLFmiT3/60zrnnHP0t3/7t608TACYdE8++aR279499rNp0yZJ0lVXXSVJuuWWW/TII49ow4YN2rp1q1555RV98IMfbOchA8CEoxYCwCjqIYCpojxZT1StVrVhwwYdPHhQfX19deds27ZNK1eudB5btmyZHn744eB+K5WKKpXK2Pq+ffsm5oABoIVOOOEEZ/2OO+7Qf/yP/1Hvfve7NTAwoPvuu0/r16/XxRdfLEm6//77tWTJEm3fvl3vfOc723HIADDhqIUAMGoy62GmXJlGY9HO42a9ah6vejnVzMRx/Ri0nWvHMi/C68wz0edCrNpsVy9a/TobEfbjwllgzJ9nz0ZhrMEIcpI0mXdugH8MrXyuhh1GMtsevvOr1ItB11kuvASNPreNDPvx4TQxQ+a9kXrnulSqjXmfG5n1RCXzuPceajALn9jn9qPLgeNV6s0LjCWlyLWG9Z4ray563/JmBs8995xmzZqlrq4uffSjH9XGjRt1+umn153b39+vefPmOY/NmzdP/f39wf2vXr1as2fPHvtZsGDBhB4/ALTa0NCQvvGNb+gjH/mIkiTRjh07NDw8rEsuuWRszmmnnaaFCxdq27Ztwf1UKhXt27fP+QGAowW1EABGUQ8BHM1a/kXjqaeeqmeffVZPPPGEPvaxj+m6667TT3/60wnb/6pVqzQwMDD28/LLL0/YvgFgMjz88MPau3evPvzhD0sa/Y8unZ2dmjNnjjOP//ACYCqjFgLAKOohgKNZy79o7Ozs1Jvf/Gade+65Wr16tc466yzdfffddef29vYWGtru2bNHvb29wf13dXWN3dX69R8AOJrcd999uvTSS3XSSSeNaz/8hxcARzNqIQCMoh4COJpNWo/G12VZ5vRUtPr6+rR582bdfPPNY49t2rQp2NMRAI52L730kr773e/qoYceGnust7dXQ0ND2rt3r/Nfrhv5Dy9dXV0tPV4AaAVqIQCMmox6mClT9pv/bw2bLoV2pOr1mhtyei96PRptv8VIL0e7ndOH0e+vqAbHIvNyZ/9qaF6xR6Pq8lv+BXs5Jn4/zNq82NVfmXlN0kaf6zBMxD7aJbc9Fb0XwumxGO3RmDY4Lzzm9nY0A5n7yiYlHT6/V2SoL2OpFJnXWM/Hwu8ceqwBLb2icdWqVfrBD36gf/u3f9Nzzz2nVatWacuWLbrmmmskSddee61WrVo1Nv+mm27SY489pi984Qv62c9+pltvvVVPPfWUbrjhhlYeJgC0zf33368TTzxRH/jAB8YeO/fcc9XR0aHNmzePPbZz507t2rWL//ACYEqiFgLAKOohgKNdS69ofPXVV3Xttddq9+7dmj17tpYuXarHH39c73vf+yRJu3btUmq+hb3gggu0fv16/fmf/7n++3//71q8eLEefvhhnXnmma08TABoiyzLdP/99+u6665TuVwrx7Nnz9b111+vlStX6vjjj1dPT49uvPFG9fX1cZdVAFMOtRAARlEPAUwFLf2i8b777ouOb9mypfDYVVddpauuuqpFRwQAR47vfve72rVrlz7ykY8Uxu666y6laarly5erUqlo2bJl+vKXv9yGowSA1qIWAsCoyaqH1TxXNc81nGeFx183bJK+w5HY81DuRjaH89pXDFkkYj1stqtGo9jhSPRI1mD8OhCJDqShfzMvNmb3H9tLJHaaxLarL/M28aPUoXkTEo+2+2x0d03Os4frHLof4zXnMPezujbOXKotJ6U0OM+JH2fuZ8O+IRLvBDunxr5xEncfeeBNlcTiyWnkeM1YYR9p4PcqnMNIyDlNmsx7t6FHIwBg1Pvf//7gH5zu7m6tWbNGa9asmeSjAoDJRS0EgFHUQwBTQcvvOg0AAAAAAABg6uOKRgAAAAAApolhVTWsXFUv9jtslquRuzg7Y961S3bM3oF62ItYVwN3nR7OG7s7tb+eRSLRobtJF+8sHR4L8eclTUWi3X2Ergbz9+1HpJt5rokW3X1SP39d2Capv+zPc94qXo48GLn2705totQ2Vp1U3VchD0WsvcN1rkj2704dy+M7EwN3ltYE3E06Mq9ehDsa647gikYAAAAAAAAA48YXjQAAAAAAAADGjS8aAQAAAAAAAIwbPRoBAAAAAJgmRvJMw7k07PWMq5pV2ytxyLs+aTivfY0w5PVetGO2L2NVXn9Fs27nZYUejUlwLNRTsdjLsf7+fPZ0NNqjMa5+T8KCQu/F2tySGWu0H2Szx+5s1+JejrHTYfsh5vZ3jPRyLBxuqX6Pwrzkn8NAL0OvN6LflzEkqVaD24TuKF/YR6gPo+T2WAz1YYzMK/RcjP3O/n4OA1c0AgAAAAAAABg3vmgEAAAAAAAAMG5EpwEAAAAAmCay3/wMe48PmyyqjTr7kWVnzLt2qZrbfdTG/H0MZzYuHYk9O88VGQvEqP31PPB4vXV3LDjkCe3D34GNs4Z3bo/Jn9doRDr+e7UwIu3tOvhUsUOIxKPzUOzZG7Nx6UJ0uGQjxoFlSYlZzzPv9QqlqvPM3Uf0F3WeLDwWikt7Me2GI+H++SiMEZ0GAAAAAAAA0CZ80QgAAAAAAABg3IhOAwAAAAAwTQzluYby3LnLtOTGnjMnAu3GK+2dpav+XaIDcelh7+7UoVh1IR4dvZt0/TF/Xuju1H4cutExqxhndkbrbuNvV9i3GbO/i3+VWCxyHXL4W7RG9AbXobHDueu0M2ZeVz8S7YyZiHXJO9tZasa819zcaNqJUcdi6TZ+7d9Z2j1Ab32cd5OOzqt3HM1F67miEQAAAAAAAMC48UUjAAAAAAAAgHHji0YAAAAAAAAA40aPRgAAAAAAponM/LiP1/qxDcn2V3S/NrA9Gwu9F02/RTuv0MvR9lR05jXXo9F2zfN7HmbBHo1+L8fwPkL8eY32TbTP5bfNa5R97mjPx8nk/P5+P8BG91H/Yf/XivVvtH0ZnbFio0szZnsZep8OO5a5Y7afo9OvsfAJM0qxvoyR6wFDfRlL7ufQ+b1ifRhjb74kUbPXJnJFIwAAAAAAAIBx44tGAAAAAAAAAONGdBoAAAAAgGmimv/mx8ub2vUsEIEe3T58vVLVXMuURebZOHM1r/+8/jw/Op0HxrLG0svRiHHe4D6K7D7DMeJGo852nv/7pw3uw99uskSfNhJ7bnhew2OBGLHciHUSmedEjv2YcrVaf568eY2+MUPxaP+4nKh3bF7k93qj6HSTuX6uaAQAAAAAAAAwbnzRCAAAAAAAAGDcWvpF49q1a7V06VL19PSop6dHfX19evTRR4Pz161bpyRJnJ/u7u5WHiIAAAAAAACACdDSHo0nn3yy7rjjDi1evFh5nuvrX/+6rrjiCj3zzDM644wz6m7T09OjnTt3jq0XMukAAAAAAGBcql4TvWZ6+VUV7qno9HyMNuILH0Osz2OjbGe8WO/FWJ/Ddgl3eZy6mnkZ8sj3RnmoX6MUvPQuT915SdbYQdnvr3L/zZYe/j4m1QQ+b0u/aLz88sud9c9+9rNau3attm/fHvyiMUkS9fb2tvKwAAAAAAAAAEywSbvrdLVa1YYNG3Tw4EH19fUF5x04cECnnHKKsizTOeeco9tvvz34paQkVSoVVSqVsfWBgQFJ0r4D2cQd/FGs+lrljSdhWhrJhtp9CEeMkXz0XBT+q9NR7vXfh3o4inqIEOrhqOlSC1/L3Jp4MK+tHzB3hLSPS9JBs91rWdUZs+uH8pGx5cHqiDNvsDo8tlwxV0cMVd1LKoZGamPD3l0qh0dq6yPOsnu81WrtbpfVkWrdZbiyivt6ZUO1/6mUV2qvUTLkXQLjjHk7rZirW+yY/ydpyLzOZl4y7H0ezXo+7P19N69tPlL7XVLvfZhUa0+QZLUDyamFY6Z6PTzwm3pY8X69irn0q2Jq4MG8sZonSYfM+20wqy1XTP2TpMpI7XMzVDU1b8T9fA2b3Y9U3eOwdc+peVX3br/umLkrtneVWp7Zq9HUFOeiMHNXaP9isTwNjyXOdmbZe648aewg7VWi/hZZZu4ubl6HzHsd7Ho2Ys6vVw+ToSQ8ZuueqY2FemjWE1Mb7bIkrx56Y6YeprY2eu+hxLxf7bL8efbvvvdvANl/V+T2TuDNvYmiVzS6b7DA403Oq/O8r/8b+XB/lyRvcfV87rnn1NfXp8HBQc2aNUvr16/XZZddVnfutm3b9MILL2jp0qUaGBjQX/3VX+kHP/iBfvKTn+jkk0+uu82tt96qT33qU638FQBMEy+//HKw1hyNfvnLX2rBggXtPgwARxlqIQCMoh4CwOHXwpZ/0Tg0NKRdu3ZpYGBA3/rWt/S1r31NW7du1emnn/6G2w4PD2vJkiW6+uqr9elPf7ruHP+KxizL9O///u+aO3duW/s77tu3TwsWLNDLL7+snp6eth3HkYBz4eJ81Bwp5yLPc+3fv18nnXSS0rSl98iaVFmWaefOnTr99NPbfo6PFEfKe+5IwLlwcT6mdi185ZVXlOe5Fi5cOK1f49fxfndxPmo4F6Omcj3k34Yu3vM1nAsX56P5Wtjy6HRnZ6fe/OY3S5LOPfdcPfnkk7r77rv1la985Q237ejo0Nlnn62f//znwTldXV3q6upyHpszZ874DnoCvX7HbXAufJyPmiPhXMyePbutz98KaZrqTW96k6Qj4xwfSTgfNZwL13Q/H1O1Fp588snat2+fJF5ji3Ph4nzUcC6mbj3k34b1cT5qOBeu6X4+mqmFk/6fZ7Isc65AjKlWq3ruuec0f/78Fh8VAAAAAAAAgPFo6RWNq1at0qWXXqqFCxdq//79Wr9+vbZs2aLHH39cknTttdfqTW96k1avXi1Juu222/TOd75Tb37zm7V37159/vOf10svvaQ/+IM/aOVhAgAAAAAAABin0q233nprq3b+zW9+U3/913+t1atX68EHH1S1WtU999yj973vfZKku+++W+VyWVdeeaUk6eGHH9bnPvc53XnnnXrooYc0d+5c/f3f/72WLl3aqkNsqVKppIsuukjl8qTd3PuIxblwcT5qOBetxzl2cT5qOBcuzsfUx2tcw7lwcT5qOBdTH6+xi/NRw7lwcT6a0/KbwQAAAAAAAACY+qbOLbQAAAAAAAAAtA1fNAIAAAAAAAAYN75oBAAAAAAAADBufNEIAAAAAAAAYNz4ohEAAAAAAADAuPFFYwusWbNGv/3bv63u7m6df/75+vGPf9zuQ2qLH/zgB7r88st10kknKUkSPfzww+0+pLZZvXq13v72t+vYY4/ViSeeqCuvvFI7d+5s92G1zdq1a7V06VL19PSop6dHfX19evTRR9t9WFPSdKxHjXzeBgcHtWLFCs2dO1ezZs3S8uXLtWfPnjYd8eS54447lCSJbr755rHHptu5+NWvfqUPfehDmjt3rmbMmKG3vvWteuqpp8bG8zzXX/zFX2j+/PmaMWOGLrnkEr3wwgttPGJMhOlYCyXqYcx0r4fUwulrOtZDamHYdK+FEvWwFfiicYI9+OCDWrlypf7yL/9STz/9tM466ywtW7ZMr776arsPbdIdPHhQZ511ltasWdPuQ2m7rVu3asWKFdq+fbs2bdqk4eFhvf/979fBgwfbfWhtcfLJJ+uOO+7Qjh079NRTT+niiy/WFVdcoZ/85CftPrQpZbrWo0Y+b7fccoseeeQRbdiwQVu3btUrr7yiD37wg2086tZ78skn9ZWvfEVLly51Hp9O5+LXv/61LrzwQnV0dOjRRx/VT3/6U33hC1/QcccdNzbnc5/7nL70pS/pnnvu0RNPPKFjjjlGy5Yt0+DgYBuPHOMxXWuhRD0Mme71kFo4fU3XekgtrG+610KJetgyOSbUO97xjnzFihVj69VqNT/ppJPy1atXt/Go2k9SvnHjxnYfxhHj1VdfzSXlW7dubfehHDGOO+64/Gtf+1q7D2NKoR6N8j9ve/fuzTs6OvINGzaMzfk//+f/5JLybdu2teswW2r//v354sWL802bNuXvfve785tuuinP8+l3Lv70T/80f9e73hUcz7Is7+3tzT//+c+PPbZ37968q6sr/+Y3vzkZh4gWoBbWUA+ph3lOLZzOqIejqIXUwtdRD1uDKxon0NDQkHbs2KFLLrlk7LE0TXXJJZdo27ZtbTwyHGkGBgYkSccff3ybj6T9qtWqHnjgAR08eFB9fX3tPpwpg3pU43/eduzYoeHhYefcnHbaaVq4cOGUPTcrVqzQBz7wAed3lqbfufjHf/xHnXfeebrqqqt04okn6uyzz9ZXv/rVsfEXX3xR/f39zvmYPXu2zj///Cl5PqYDaqGLekg9lKiF0xX1sIZaSC18HfWwNfiicQL93//7f1WtVjVv3jzn8Xnz5qm/v79NR4UjTZZluvnmm3XhhRfqzDPPbPfhtM1zzz2nWbNmqaurSx/96Ee1ceNGnX766e0+rCmDejSq3uetv79fnZ2dmjNnjjN3qp6bBx54QE8//bRWr15dGJtu5+IXv/iF1q5dq8WLF+vxxx/Xxz72Mf3Jn/yJvv71r0vS2O883T83Uwm1sIZ6SD18HbVweqIejqIWUgst6mFrlNt9AMB0s2LFCj3//PP64Q9/2O5DaatTTz1Vzz77rAYGBvStb31L1113nbZu3cqXjZhQ0/3z9vLLL+umm27Spk2b1N3d3e7Dabssy3Teeefp9ttvlySdffbZev7553XPPffouuuua/PRAa1FPaQevo5aiOmMWto11BQAAAPoSURBVEgttKiHrcEVjRPot37rt1QqlQp3ZNqzZ496e3vbdFQ4ktxwww369re/re9///s6+eST2304bdXZ2ak3v/nNOvfcc7V69WqdddZZuvvuu9t9WFMG9Sj8eevt7dXQ0JD27t3rzJ+K52bHjh169dVXdc4556hcLqtcLmvr1q360pe+pHK5rHnz5k2bcyFJ8+fPL/zHjCVLlmjXrl2SNPY7T+fPzVRDLRxFPaQeWtTC6Yl6SC2UqIU+6mFr8EXjBOrs7NS5556rzZs3jz2WZZk2b95M77lpLs9z3XDDDdq4caO+973vadGiRe0+pCNOlmWqVCrtPowpYzrXozf6vJ177rnq6Ohwzs3OnTu1a9euKXdu3vve9+q5557Ts88+O/Zz3nnn6Zprrhlbni7nQpIuvPBC7dy503nsX/7lX3TKKadIkhYtWqTe3l7nfOzbt09PPPHElDwf08F0roUS9dCiHtZQC6en6VwPqYU11EIX9bBF2nsvmqnngQceyLu6uvJ169blP/3pT/M/+qM/yufMmZP39/e3+9Am3f79+/Nnnnkmf+aZZ3JJ+V//9V/nzzzzTP7SSy+1+9Am3cc+9rF89uzZ+ZYtW/Ldu3eP/bz22mvtPrS2+LM/+7N869at+Ysvvpj/8z//c/5nf/ZneZIk+Xe+8512H9qUMl3rUSOft49+9KP5woUL8+9973v5U089lff19eV9fX1tPOrJY+8smOfT61z8+Mc/zsvlcv7Zz342f+GFF/K///u/z2fOnJl/4xvfGJtzxx135HPmzMn/9//+3/k///M/51dccUW+aNGi/NChQ208cozHdK2FeU49fCPTtR5SC6ev6VoPqYVx07UW5jn1sFX4orEF/uZv/iZfuHBh3tnZmb/jHe/It2/f3u5Daovvf//7uaTCz3XXXdfuQ5t09c6DpPz+++9v96G1xUc+8pH8lFNOyTs7O/MTTjghf+9738uXjC0yHetRI5+3Q4cO5X/8x3+cH3fccfnMmTPz3/u938t3797dvoOeRP4/JqfbuXjkkUfyM888M+/q6spPO+20/N5773XGsyzLP/nJT+bz5s3Lu7q68ve+9735zp0723S0mCjTsRbmOfXwjUznekgtnL6mYz2kFsZN51qY59TDVkjyPM8n48pJAAAAAAAAAFMXPRoBAAAAAAAAjBtfNAIAAAAAAAAYN75oBAAAAAAAADBufNEIAAAAAAAAYNz4ohEAAAAAAADAuPFFIwAAAAAAAIBx44tGAAAAAAAAAOPGF40AAAAAAAAAxo0vGgEAAAAAAACMG180AgAAAAAAABg3vmgEAAAAAAAAMG7/H6PpXzyACRhWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1600x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Testing\n",
    "scale = 20;\n",
    "img_sc_nn = nnInterpolation(image, scale);\n",
    "img_sc_bl = bilinearInterpolation(image, scale);\n",
    "img_sc_bc = bicubicInterpolation(image, scale);\n",
    "\n",
    "\n",
    "fig=plt.figure(figsize=(16, 4), dpi= 100, facecolor='w', edgecolor='k')\n",
    "plt.subplot(141)\n",
    "plt.imshow(image)\n",
    "plt.title(\"Original image\");\n",
    "\n",
    "plt.subplot(142)\n",
    "plt.imshow(img_sc_nn)\n",
    "plt.title(\"Nearest neigh\");\n",
    "\n",
    "plt.subplot(143)\n",
    "plt.imshow(img_sc_bl)\n",
    "plt.title(\"Bilinear\");\n",
    "\n",
    "plt.subplot(144)\n",
    "plt.imshow(img_sc_bc)\n",
    "plt.title(\"Bicubic\");\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
